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TECHNICAL PAPER 


NONLINEAR OPTIMIZATION WITH LINEAR CONSTRAINTS 
USING A PROJECTION METHQD 

L INTRODUCTION 


The solution to the extremization problem with a nonlinear objective and Unear equality and inequal- 
ity constraints has application to many fields of science and business. Extremization means either maximiza- 
tion or minimization of an objective function with respect to a set of decision variables that may be required 
to satisfy a set of equaUty and inequality constraints. Finding a solution that simultaneously satisfies the 
constraint equations and extremizes the objective function is, in general, not a straightforward procedure. 
To this end, many methods of extremizing functions have evolved. 

Problems can be divided into two broad categories, linear problems and nonlinear problems. The 
Unear problems have a Unear objective function and linear constraint equations which can be expressed in 
the following general form: 


extremize f(xj,X 2 ,...,Xj^) = CjXj + C 2 X 2 + ... + 


subject to 


n 

^aijXj-bi<0 , k 

j-1 


n 

L ^ij^j - bi = 0 , i = k+1, k+2, h 
j=l 


n 

^ ajjXj-bj>0 , i = h+l,..., m 

j=l 


If the variables appear nonlinearly in either the constraints or the objective function, then the problem is 
considered to be nonlinear and, consequently, is generaUy more difficult to solve. Within this class of prob- 
lems consider these two subclasses: those where the nonUnearities appear only in the objective function, and 
those where the nonlinearities can appear in both the objective and the constraint equations. 

This paper addresses the problem described by a nonlinear objective function and linear constraints. 
The solution technique is based on a method proposed by Rosen (1960) [1] called the Gradient Projection 



Method. In this method, if any of the constraint equations are violated during the unidirectional search, a 
projection method is used to generate a new feasible direction of search. The functions considered here will 

be convex within the region of interest. The functions wiU also be of class (first and second order partials 
exist and are continuous). The feasible region is defined by a convex polyhedron formed by the linear con- 
straint equations. Rosen’s method, unfortunately, requires the inverse of the matrix formed from the 
normals cf the binding constraints. Binding constraints are the equality constraints and the violated inequal- 
ity constraints that are treated as equality constraints. Rosen attempted to reduce the complexity of this 
problem by updating the inverse matrix with a recursive method rather than recalculating it whenever there 

T 1 

was a change in the binding constraint set. This recursive method depends on the knowledge of (N*N)'^, 
where N is the matrix of the normals of the binding constraints and is the matrix transpose [ 1 ] . Only 

one hyperplane can be subtracted or added to the projection matrix^ at a time. Table 1 shows the number 
of computations needed to update the projection matrix for the subtraction of one constraint by the Rosen’s 
method and the method proposed in this paper. This table does not include the computations that would be 

required by the calculation of (N N)*^. The dimension of the projection matrix calculated by either method 
is equal to the number of independent variables. The Rosen method does not calculate the projection matrix 
with' matrices of this dimension. The rank of the matrices of this method varies with the number of 
independent binding constraints. The proposed method saves operations by taking advantage of the matrix 
symmetry. The last column in Table 1 shows the ratio of the number of operations required by the two 
methods. This paper proposes a solution technique that does not require the calculation of the inverse of 
the matrix of the normals of the binding constraints and, for a large number of problems, requires a smaller 
number of computer operations than does Rosen’s method. The proposed method provides similar results 
with fewer computations, thus in general reducing the computer time. This method is able to add or remove 
many constraints at a time which the Rosen method does not. A comparison between the two methods is 
seen in Figure 1 for the case where there are 20 independent variables. 

Section II of this paper discusses the general nonUnear optimization problem subject to linear con- 
straints. In Section III the problem of finding an initial feasible point is addressed. Phase I of the Two-Phase 
Simplex Method is used to provide a feasible point when a user-provided starting point is infeasible or when 
the problem includes equality constraints, since these must always be satisfied. Section IV discusses a tech- 
nique for locating an extrema within unconstrained feasible space. The method is the Davidon-Fletcher- 
Powell (DFP) which is a variable metric technique [ 2 ] . Section V analyzes the solution to the problem with 
constrained extremum and presents the main thrust of the paper. The application of the gradient projection 
method to locating constrained extrema wiU be discussed. The differences between the method proposed in 
this paper and Rosen’s method will also be detailed. Section VI contains the conclusions. Appendices A and 
B provide the user’s guide and the results of the test cases to which the program was appUed. Appendix C 
contains a description of the program used in this paper and Appendix D contains the program listing. 


II. STATEMENT OF THE PROBLEM 


Let Xj,i=l,2,...,n, be the coordinates of the point x in n-dimensional Euclidean space. Points in space 

wiU be defined with superscripts, and elements of points will be defined with subscripts. Vectors wiU be 
represented by columns of elements as 


T 

1. An n X n real matrix P is caUed a projection matrix if and only if P = P and PP = P. 


2 



TABLE 1. COMPUTATION COMPARISON FOR UPDATING PROJECTION MATRIX 



2 4 6 8 10 12 14 16 18 


BINDING CONSTRAINTS 

Figure 1 . Operations to update the projection matrix. 
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^2 


T 

X will be defined as the transpose of the x vector. In this notation, the general maximizing nonlinear pro- 
gramming problem with hnear constraints can be expressed as: 


max: f(x) = f(xj,X 2 , Xj^), 


subject to the hnear equalities and inequalities of the form: 


(a^)^ X - b| = 0 , i = 1 ,2,... ,k<n 
(a^)^x-bi>0 , i = k+l,...,p . 


The vectors a\ i=l,2,...,p, provide the constraint equation coefficients and the bj, i=l,2,...,p are scalars. The 
symbol g^ will represent the ith constraint. These constraints restrict the solution to k hyperplanes and p-k 

half spaces. Their intersection, R, is a convex polyhedron called the feasible region. R consists of all the 
points that lie on the equaUties and within the half spaces. Points in the equality constraints will lie on the 
boundary of this region. The constraints are assumed to be linearly independent. The problem is depicted in 
Figure 2. The curves hues in Figure 2 are the level contours of f(x). This figure shows the linear constraints 
S 2 ’S 2 ’S 4 < 0 and gj,g 5 > 0 which outline a feasible space that contains x*. Starting at the unconstrained 

point x°, the search follows the gradient of f(x°) until a constraint is encountered at x^. The projection 

search vectors p^ and p^ are calculated at x^ and respectively, as new constraints are encountered at 
these points. A constrained optimum may be found at points that are not located at vertices in nonlinear 
problems. Inspection of this figure shows that the optimum point x* occurs on constraint gj, away from a 
constraint vertex. 


III. FINDING A FEASIBLE POINT 

Most iterative search methods require an initial starting point. These points must either be supplied 
by the user or be generated by an initial point algorithm in the method. For linearly constrained problems 
the Simplex method can be used to find an initial feasible point if one exists. Many nonlinear methods are 
developed based on a quadratic function and extended to the general nonquadratic case by approximations 
of the quadratic. For initial points which are not in the vicinity of the optimal solution, the function may 
not be adequately represented by a quadratic approximation and the solution algorithm generally takes 
longer to converge. If the user supplied initial point for linear constrained problems is not in R, then one 
method of obtaining an initial feasible point is to use Phase I of the Simplex method. The initial point will 
then lie on the boundary of R at one of the vertices formed by a subset of the constraints, and will auto- 
matically satisfy the equality constraints if they are consistent. 
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The Simplex method extremizes a hnear objective function subject to the constraint set 


Ax = b 


( 1 ) 


and 


Xi> 


i= 1,2,. ..,n 


where A is an m x n matrix and b is an m-dimensional vector. 

The null space of an n x s matrix C is the subspace of all s-dimension vectors y such that 


Cy = 0 


Denote the subspace of the vectors that have this property by N(C). If x and y are members of this set, then 
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C(x+y) = Cx + Cy = 0 , 


and if k is a scalar then 


C(kx) = k(Cx) = 0 . 


Let N(C) be the null space of C and have dimension q. Then the set has q linearly independent 
vectors which form the basis for the set defined by C. The column rank of a matrix is equal to the number 
of columns minus the null space dimension. Let the basis of the null space of C be 

v^,v^,...,v^ . 


This set can be extended to a basis for an s-dimensional space by adjoining wl, j=l,2.,„.r vectors that are 
hnearly independent 


v^,v^,...,v‘i,w^,...,w^ 


where 


q+r = s . 


Now every column vector x in C can be written 


q r 

x = ^ a^vi + ^bjwJ 
i-1 j=l 

Therefore, 


q r q r 

Cx = C ^ aj vi + C ^ bj wj = ^ aj (Cvb bj (Cwi) 
i=l J=1 i=l j=l 


Since Cv^ = 0 then 
r 


Cx = ^bj(Cwi) 

j=l 
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Hence, the r vectors Cw^ span the co-domain. These vectors are also linearly independent. If 


ki(Cwb + ... kj.(Cw*’)= 0 , 

then 

Ch=0 , 

where 

h = kjw^ + ... + kj-w*- . (2) 

So h belongs to N(C) and is a linearly independent combination of the r vectors 
q r 

h=^l/ + X;kjwj , 

i=l j=l 

a contradiction. 

The Simplex method of finding an optimal solution is an iterative update process that solves the 
constraint equations for their vertices. The direction of the move from one vertex to the next is done in a 
manner that increases or decreases the objective function depending on whether it is a minimization or 
maximization problem. Letting z^ be the value of the objective function for the previous iteration, 

z = Zq + f (optimizing vector) . 

The optimizing vector transforms the last value of the objective function z^ to a new value z that is closer 
to the optimal value. Let A be partitioned into matrices B and N such that 

A = (B,N) 

where B is an m x m invertible matrix and N is an m x (n-m) matrix and the rank of A(B,N) is m. Let 
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where the partitions of x have the dimensions m and n-m, respectively, and where Xg is associated with the 
B partition of A and xj^j with the N partition of A. From (1) 


Ax = (B,N)x = Bxg + Nxj^ . 


A basic solution will be given by x if 


Xg = B‘^b 
Xn= 0 . 


If X > 0 , then this represents a solution that satisfies all constraints and is called a basic feasible solution 
(BFS). In general, only a small portion of the vertices will have to be examined in order to locate the optimal 
solution to a linear objective function using the Simplex algorithm. The Simplex method requires an initial 
feasible point. Given this, each succeeding vertex will also be feasible. When a better objective function 
value cannot be found at any adjacent vertex and all xgjj > 0, i=l,2,...,m, an optimum has been achieved. 
Consider the following minimization problem : 

T 

mm z = c X 


such that 


Ax = b 


x>0 . 


Given a basic feasible solution Xg^, 


~^b’ 


"B'^bl 

_Xn_ 


1 

o 

1 


let cg be the coefficients in the objective function associated with the xg vectors and cjq with the xj»j 
vectors. The objective function z can be given by 

B-^b 

= c^B‘lb . 

0 


T T T 
Z = C X = (Cg,cj|) 
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Let 

X^= (Xb,Xj^)^ 

be an arbitrary basic feasible solution. Then 
Xg > 0 , = 0 , 

and 

Ax = Bxg + Nxj^ = b 
Premultiplying by and rearranging 

Xg = B'^b - B'^Nxjvj 

or 

xg = B’^ b- ^ B'^a^Xj (3) 

jeD 

where aP is the jth column vector in A and where D is the current set of indices of the nonbasic variables. 
The objective function is 

z = Jx = Cg^xg + = c^(B-^b - Y, B'^a^Xj) + Y ^jXj 

jeD jeb 

Let Zq represent the objective function at some arbitrary basic feasible solution. Then 

z = Zq - E (^j • ‘^j)^j ’ 

jeD 
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where Zj = cg^B’^ai for each nonbasic variable. Since z is to be minimized, it can be seen from (4) that 
whenever ^ - cj > 0 the objective function decreases by introducing some nonbasic variable into the basis 

at a positive value. This rate of decrease will be the largest if we pick the most positive of these zj - cj, say 
Zk - cjj. From (3) 

xg = B' ^ b-B‘ ^ a^xj^ = b-y^xj^ , 


where b = B'^b and y^ = B'^a^. Then as xj^^ increases positively from zero, the basis is modified as shown in 
vector form 


’^Bl ' 


’^1 ' 



^B2 

= 

^2 

- 


_^Bm_ 


1 

s 

1 


.Vm’l 


Xk . 


(5) 


If < 0 i=l,2,...,m, then xgj increases as Xj^ increases indefinitely. If all y^^ > 0 then xj^ can increase 
only until one of the xgj = 0 due to the x > 0 constraints. In the absence of degeneracy (degenerate solu- 
tions are those where the value of at least one of the basic variables equals zero), bj. > 0, r=l,2,...,m and then 
Xk = bj./yj.^ > 0 . From (4) 


2=^o-(^k-CkH • 

From the choice of Zj^ - cj^ > 0 for a minimization problem. 


( 6 ) 


Z<Zq . 

Substituting (6) into (5) 

xb = b. - (yi^/yr^)br , i=l,2,...,m . 

All other xj = 0 . From above, xgj. = 0 becomes a nonbasic variable, and Xj^ = bj./yj.^ has been added to 

replace it, forming a new BFS. This new BFS will decrease the value of the objective function and will also 
satisfy the constraints. 

Phase I of the Simplex method requires the problem to be set up in the standard form [3] . This is 
accomplished by the use of slack and artificial variables. Standard form is obtained when all variables are 
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non-negative, ail of the constraints are equations, and all constants on the right-hand side are non-negative. 
In this form, row column operations can be performed. In some problems, the coefficient matrix A may not 
contain a full unit matrix of the same order as the number of constraints. This occurs when there are 
equality and less than or equal to constraints. In this case no BPS can be obtained. Artificial variables can 
be added to the constraints to augment A such that a BPS exists. Non-negative slack variables are intro- 
duced as necessary to transform the set of less than or equal to constraints into equalities. This set is then 
examined to determine if a full unit matrix is present. If it is not, non-negative artificial variables are added 
to the appropriate equations to create a fuU unit matrix. The Simplex method begins at the origin of the 
primary decision variables with all of the vector b allocated to the artificial and slack variables. These arti- 
ficial variables must be driven to a zero value in the final solution. In Phase I of the two-phase Simplex 
method, the objective function is replaced with an objective function that sums the artificial variables. This 
objective function will be minimized. Consider the following set of constraints: 


Xj - 2x2 + X3 < 11 


-4xj X2 + X3>3 


2 xj - X 3 = - 1 
xj>0 i= 1,2,3 . 


The addition of non-negative slack variables X 4 and xj allows the problem to be expressed in standard form: 


Xj - 2x2 + X3 + X4 = 11 

-4xj + X 2 + 2 x 3 - X 3 = 3 

- 2 xj + X 3 = 1 

xj> 0 i= 1,2,3 ,4,5 . 


Artificial variables can be added to the last two equations to produce the required unit matrix. The Phase I 
Simplex method will reduce these artificials to a value of zero, if a BPS to the problem exists, and thus 
remove them from the solution. Prom the previous example, after adding the artificial variables xg and x^, 
the Phase I problem is 


Min z = xg + x -7 , 


subject to 

xj - 2x2 + X3 + X4 = 11 
-4X1 +X2 + 2X3-X5 + Xg = 3 


I 


11 



-2xj + X3 + Xy = 1 


Xi>0 i= 1,2, 3, 4, 5, 6, 7 . 


This problem can be represented in the following tableau form: 




0 

0 

0 

0 

0 

1 

1 



Basis 


^2 

^3 

X4 

^5 

^6 

^7 

b 

0 

X4 

1 

-2 

1 

1 

0 

0 

0 

11 

1 

^6 

-4 

1 

2 

0 

-1 

1 

0 

3 

1 

X? 

-2 

0 

1 

0 

0 

0 

1 

1 

c row 

6 

1 

-1 

-3 

0 

1 

1 

0 

0 



The tableau is a representation of the problem in detached coefficient form. The Cg column is made up of 

the objective function coefficients of the present basic variables listed in the basis column. The coefficients 
of the objective function are shown in the Cj row. The columns contain all of the coefficients associated with 

the variables with the exception of the c row. The columns that contain the present basic variables also con- 
tain the unit matrix. The bottom row is called the cost row and ^ = Cj - (cgaj). The b column represents 

the value of the basic variables listed on the same row. The program presented in the last section of the paper 
addresses the problem of finding an initial feasible point. The most efficient way to accomplish this is 
through the user’s organization of the constraints and variables. Nonstandard forms can usually be put into 
standard form by simple substitutions or transformations. For example, in the case where 


-00 < Xj < 00 , 


the substitution of 



where 


.+ x-~ 
1 ’ ^1 


> 0 


can be made. The scalar expression 
|ax| < b 
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can be expanded in the form 


-b < ax < b 


which can be used in the form 


ax < b 


-ax < b . 


Many other nonstandard forms can be standardized with similar substitutions or transformations. 


IV. UNCONSTRAINED OPTIMIZATION 


Section III discussed the problem of finding an initial feasible point. An initial feasible point is 
essential to using the method proposed by this paper. Once a feasible point is found, a nonlinear optimiza- 
tion method is then applied to the problem. The method proposed here uses the DFP method when the 
present point is not located in any constraint hyperplanes and does not violate any of the constraints. 

The advantage of combining the Rosen gradient projection method (1960) [ 1 ] with the DFP method 
first introduced by Davidon (1959) and later modified by Fletcher and PoweU (1963) [2] is seen in the rate 
of convergence to the optimal solution. The DFP method, originally developed as a nonlinear unconstrained 
technique [4,5,6] , was adapted to a linear constraint method by Goldfarb (1969) [7] . The DFP method is 
considered to be one of the most powerful of the nonlinear search techniques and was developed to solve 
the quadratic function. The method becomes iterative when extended to nonquadratic nonlinear functions. 
By using the Taylor series, it is possible to determine the value of a function in the neighborhood of a known 
point. As the neighborhood becomes small, the Taylor series approximates a quadratic function. The smaller 
the neighborhood the better the approximation. The DFP method combines the best features of the steepest 
descent technique with those of the Newton method and with few of their drawbacks [8]. The steepest 
descent technique generates a search vector using knowledge of the first partials of the objective function. 
The negative gradient of the objective function is represented by the vector 


-Vf(x) = - 


9f(x) ~ 

9xj 


9f(x) 


3x 


n J 


(7) 


and points in the direction of the most rapid decrease of the objective function from the point x. Steepest 
descent converges to the optimum solution slowly when the function is highly nonlinear. A Newton method 
converges in a single step when operating on a quadratic objective function. This method requires that the 
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matrix of second partials be calculated and inverted at each iteration. The DFP method approximates the 
Hessian inverse using only first partials, and improves this approximation at each step. This method gener- 
ates the inverse Hessian after n steps in an n-dimensional quadratic problem. Like the method of conjugate 
gradients, the algorithm of the variable metric DFP method is designed to extremize the foUowing function 

about an arbitrary point x°: 


f= l/2(x-x°)^A(x-x°) 


( 8 ) 


by conducting a sequence of one-dimensional line searches that begin at an arbitrary point Xj and locate 
improved points by the relation 

xi+1 = x^ + a^p^ , (9) 

where oj is a positive scalar constant and p^ is a search vector. The constant oj is chosen so that is 

located at the extremum along the direction of search. Since this extremum is located at a point x^"^^ where 
the gradient of the objective function is perpendicular to the search vector, their dot product is zero [8]. 
Let 


y(a|) = f(x^ + a-pi) 


9y(ai)/3ai = (p^)^Vf(x^'^^) (10) 

= 0 . 

The extremum along the direction of search occurs at the point where the search vector is tangent to 
a contour of the objective function as shown in Figure 3. From (8), 

Vf(x^) = A(x^ - x°) . (11) 
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Using (9) recursively yields 


n-1 

x^ = x^+J2 OjP^’ • (12) 

j=l 


Subtracting x° from (12) and premultiplying by A, 


n-1 

A(x”-x°) = A(xl-x°) + OjApj . (13) 

j=l 


Then from (1 1) and (13) 

n-1 

Vf(x”) = Vf(xl) + ^ OjAp^ 

j=l 

or as a special case 

Vf(xl"'’l) = Vf(xl) + OjApl . (14) 

If a p^ can be chosen so that it is A conjugate to all other search vectors, i.e., 


(p^)’^Apj = 0 , i¥=j , 


(15) 


then two important results can be concluded. First, all of the search vectors p will be linearly independent, 
and second, the quadratic form wiU be extremized in no more than n steps in n-dimensional space [9]. 
Define 


pi = -HVf(xl) , i = 0,l,...,n-l , 


(16) 


where H is an n x n symmetric positive definite matrix. Fletcher and Powell [2] proved that the iterative 
search method developed by Davidon satisfies the condition specified in (15). 

This method will be developed and some of its interesting properties discussed. H will always be 
symmetric and positive definite if Hq is chosen to be symmetric and positive definite and if is calcu- 
lated by 

15 


■ 


+ q , 


(17) 


where and Cj are n x n symmetric matrices that are calculated at each step. Using (17) recursively, 


n-1 n-1 

Hn = Ho+E Bj+E Cj 

j=0 j=0 


(18) 


By choosing 


n-1 


j=0 


and 


n-1 

E Cj=-H„ , 
j=0 

then (18) reduces to is usually chosen to be I to assure that it is positive definite and 

symmetric. The construction of B shows the importance of the A conjugacy as expressed in (15). The 
individual Bj’s in (18) can be determined by this A conjugate condition. Let T be an n x n matrix with 

columns consisting of the search vectors p^, i=0,l ,...,n-l . Then 

T^AT=D , (19) 


where D is an n x n diagonal matrix. This follows from (15) since 

dy = (P^)^Api = 0 , i^j . 

The double subscripting of a variable will be taken to denote a matrix element. Since A is positive definite, 
A'^ exists and A can be solved from (19) 

A = (tT)-1dT"* = (TD-^T^)-! . 
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Hence 


A"1 = TD-It"^ 


and thus. 


n-1 

Bj = TD'1tT . 
j=0 


From linear algebra, the diagonal terms of D'^ are 1/djj. Therefore, 


n-1 n-1 

Z Z (l/djj)p’(p’)^ (20) 

j=0 j=0 


where 


djj = (p^)^Apj 


Substituting this last expression into (20) gives 


Bj = (pj(p^)^)/((pj)^Apj) 


( 21 ) 


and substituting (14) into the denominator. 


n-1 n-1 

E Bj=E aj(Ap’)h/((p’)'^(Vf(xj+l)-Vf(xj)) 
j=0 j=0 


or, in corresponding terms, let 


Bi = ai(p^pb’^)/((pb'^( Vf(xi+1) - Vf(xi))) 


( 22 ) 


17 



for i=0,l . In a similar manner, will be developed. Post-multiplying (17) with Ap\ we have 


Hi+iApi = HiApi + BiApi + qApi . 

Substituting (21) into the middle right-hand term, 

q+1 Ap^ = HjAp^ + (p\p^)^Ap^)/((p*)’^Ap^) + qAp^ 

= H^Ap^ + p^ + qAp^ . (23) 

If q is chosen so that H^+jAp^ = p\ then p^ is an eigenvector of Hj.j.jA. Fletcher and Powell used this 
property to show quadratic convergence [2] . Using this result in (23) 

qAp^ = -HjApi . 

Since (p^)^A^HjAp^ is a scalar value, then 

1 = ((pi)TATHiApi)/((pb'^ATHiApi) . 

Hence, 

qApi = -qApi((pi)TATHiApi)/((p')'^ATHiApi) 

or 

q = -(qApi(pi)TATHiT)/((pi)TATH.Api) . 

The last term in the numerator is a result of the symmetry of H. Again, using (14) 

q = -(q(Vf(xi‘‘’l) - Vf(xi))( Vf(xi‘*‘l) - Vf(xi))TqT)/((^f(xi+l) - Vf(x^))Tq( Vf(xi''’l) - Vf(xi))) . 
The third term in the numerator follows from the symmetry of A and from the linear algebra identity 



The search algorithm of (9) with the a search vector defined by the relationship 


pi = -HiVf(xi) , 


where Hj is updated according to (18) and where the exact step Oj is found along each search vector, will 

converge to the extremum of a quadratic in no more than n steps in n-dimensional space. Examination of 
the Taylor’s series expansion of a nonquadratic nonlinear function around the optimal solution shows that 
as the solution approaches the optimal solution, the function becomes dominated by the quadratic terms. 
In the program presented in this paper, the constant oj is only an approximation of the exact distance from 

x^ to the extremum along p^. The usual method for calculating (Xj is to curve-fit the points that bracket or 

lie in the vicinity of this extremum. This program uses a quadratic fit. No effort was made to determine the 
consequences of this error on the rate of convergence. The H matrix is reset to I, since H may tend to 
accumulate errors due to numerical truncation and also due to the errors caused by not finding the exact 
extremum point along each search vector. In this latter error source, the shght inaccuracy causes the dot 
product between the search vector and the gradient of the function at that point to differ from zero, and it 
is also generally necessary to use double precision to reduce the truncation problem. These errors will 
accumulate in the H matrix at each update. A decision must be made to reinitialize the H matrix after a 
number of updates. The usual practice is to reset the H matrix to I after n steps in an n-dimension space. 
The first search vector after this update is the negative gradient of the objective function. Hence, after every 
n steps when H = I, there is a search taken in the direction of steepest descent. 

The advantages of the DFP method lie in the convergence properties that approximate the Newton 
method while using only first order information similar to the steepest descent method. H only approxi- 
mates closely after n steps. The combining of the unconstrained search with the projection technique by 
Rosen allows a search of the feasible region R. Unlike linear programming methods, which are constrained to 
examine only the vertices, this method is free to traverse interior feasible space. 


V. CONSTRAINED OPTIMIZATION 


Section IV discussed a method of nonlinear optimization for problems that have no constraints. 
Many nonlinear problems do have linear constraints, however, and this section will address this type of 
problem. 

The projection method is used at the boundaries of a convex polyhedron where the extremum lies 
outside or on the boundaries of the feasible region R. This method is based on the projection of the gradient 

of a nonlinear function onto these boundaries [ 10] . The space E*” is defined as a Euclidean space of dimen- 
sion m and contains a subspace that is spanned by the normals of q binding constraints. The normals of these 
binding constraints are the column vectors of a matrix Ng, where the subscript q denotes that the matrix is 

an m X q matrix. The symbol gj will designate the ith hyperplane associated with the ith constraint. The 
matrix 
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II HIM II II II nil 


is used as pseudo-inverse of the matrix Nq. Let a set of q linearly independent unit vectors u^, i=l,„.,q < m, 
be normals to a set of q linearly independent hyperplanes. Let 

Nq = (u^,u^, . (24) 

Because of the linear independence of the q vectors, the q x q symmetric matrix Nq^Nq is nonsingular, and 
therefore its inverse (Nq^Nq)"^ exists. Let Q be a subspace of E*’’ that is spanned by the vectors u^,u^,..., 

u^l. Let V be the orthogonal complement of Q in E*”. Then V is a subspace of dimension (m-q). Rosen’s 
projection method uses a matrix given by 

fq = Nq(N jN,r>NqT (25) 

that projects a general vector in E™ onto Q. To show this, let v e V. Since V is orthogonal to Q, 

(u^)^v = 0 , i=l,..., q 


or 


NqV = 0 . 

Using this in (25), 

PqV = 0 . 

In Q let w be a vector, which can then be represented as a linear combination of the q vectors that span Q. 
Then 


q 

w=23aigi , 

i=l 


(26) 


and from (24), 


w = Nqa , 
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where the vector a is the column vector of the coefficients in (26). Therefore, 


= N,a 

= w . 

These results show that Pq is a projection matrix that takes any vector in E™ into Q. Let Pq be defined by 



If V and w are defined as before, then 
PqV=v-PqV=v-0=v 

and 

PqW = w-PqW = w- w = 0 . 

From this, any arbitrary vector, x e E™, decomposed as 
X = V + w 

is projected onto the intersection of the q hyperplanes. It is also clear that if q = m, then Pq “ I 
PqX = 0 . 

From vector analysis, the projection of a vector z onto the ith hyperplane where V nj(x) is the normal to 
this hyperplane is 

p^ = z - d^ Vn^(x) . 
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The variable is a scalar constant. Figure 4 depicts this where z = V f(x). If p is specified to be the projec- 
tion onto Q, then 


q 

p = Vf(x) - X) <ii ^i^i • 
i=l 


The projection p will lie along the intersection of the q hyperplanes. 

A Langrangian function L is a scalar function of x and X where, if X^ is a scalar multiplier and h|(x) 
designates the ith equality constraint, then 


5 e 

L(x,X) = Vf(x) + X) \ Vn^(x) + X ^i • 

i=l i=s+l 


At a constrained extremum, the gradient of the Lagrangian vanishes [8,1 1 ] and since Xj = 0 for nonbinding 
constraints. 


s C 

VL = Vf(x) + X \ ^nj(x) + X \ Vhj(x) = 0 . 
i=l i=s+l 


Vf(x) 



22 


r 


Upon substituting the d^’s for X;^’s, the equivalency between the Lagrangian and the projection equation is 
obvious. The Kuhn-Tucker [12] conditions show that, if the d^’s have the proper signs, then 

s £ 

p = Vf(x) - XI di Vnj^(x) - Y. ^hi(x) = 0 

i=l i=s+l 

at a bounded optimum. Thus, the projection vector p vanishes at this point. The Kuhn-Tucker conditions are 
based on Farkas Theorem [12,13] which states that only one of the two following systems has a solution: 

System : 1 Ax < 0 and cx > 0 

or 

System : 2 w^A = c and w > 0 

where A is a given m x n matrix, c is a given vector of dimension n, and x,w are variable vectors of dimen- 
sion m and n. The first part of the proof of this will be by contradiction. If System 2 has a solution w such 

that w"*"a = c and w > 0, let x be such that Ax < 0. Then c"*^x = w"*^Ax < 0 which contradicts the given 
statement. If System 2 has no solution, then c ^ S = (w"^A:w>0). There is some x such that c'^x > w"^Ax 
for all w > 0. If w = 0 then c"^x > 0, and as w approaches infinity Ax < 0 completing the proof. The Kuhn- 
Tucker conditions follow from this theorem [14]. Let Z be a nonempty open set in and let f:E™^E^, 
n^:E^’^E^ for i=l,...,s, and h|:E™->-E^ for i=s+l,...,l. Consider the problem L where n^(x) is the ith 
constraint equation. 

Min f(s) 
such that 

ni(x)<0 , i=l,2, ..., s 

hi(x) = 0 , i = s-l-1, ..., 1 

and 

X e Z . 

Let X be a feasible solution, and let I = (i:nj(x) = 0). Suppose that f and n- for i ^ I are differentiable at x, 
that n^ for i e I is continuous at x, and that h^ for i=s+l,...,l is continuously differentiable at x. Further 

suppose that Vn.(x) for i e I and Vhj^(x) for i=s+l,...,l are linearly independent. If x* solves problem L 
locally, then there exist scalars for i= !,...,! such that 


23 



and 


Vf(x*)+^ Xi*Vni(x*)+ Xi*Vh|(x*) = 0 
iel i=s+l 


\*>0 , iel . 

In addition to the above assumptions, if n^ for i y I is also differentiable at x, then the Kuhn-Tucker condi- 
tions could be written in the following equivalent form : 


s S. 

Vf(x*) + 2^ Xi*Vn^(x*)+ ^ Xi*Vhi(x*) = 0 
i=l i=s+l 

\*nj(x*) = 0 , i=l,...,s 

Xi*>0 , i=l,...,s . 


This shows that Vf(x*) is a hnear combination of the normals to the constraints. These normals describe a 
cone which contains the function gradient. If A is the matrix of the constraints, then by Farkas Theorem, 
the intersection of the half space defined by Ax < 0 and the open half space cx > 0 is empty and therefore 
has no solution. A geometric interpretation is shown in Figure 5. If some d^ = 0 or has the wrong sign, then 

the associated constraint hyperplane is nonbinding and must be dropped from the set used to define the 
projection matrix. At x* the gradient of the Lagrangian vanishes when Vf(x*) becomes a linear combination 
of the binding constraints. The projection scalar d and the Lagrange multiplier X are equivalent. 



Figure 5. Graphical display of Farkas Theorem. 
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The number of computations to determine (Nq^N^)"^ is formidable when is very large. If this 
inverse was recalculated each time a change was made in the defining set of g^’s, the required computer time 
could become prohibitive. Rosen used a recursive scheme to update the (Nq^Nq)“ matrix rather than 

recalculating it each time a change was made in the binding constraint set. Although this method greatly 
reduced the amount of computation necessary to calculate a new inverse, it can remove or add only one 
constraint at a time. 

In this paper, a method is presented that projects an arbitrary vector x c onto a subspace Q C 
without the necessity of calculating an inverse. For a projection onto a hyperplane, it is necessary to find 
an exact point in the constraint hyperplane at which to calculate the gradient of the objective function. 
This gradient will then be projected onto the intersection of the binding constraints. The development of 
this method follows. 

A search along the vector Vf(x^) is made using a scalar multiplier Let 

xi+l = xi + oirf Vf(xi) (27) 


where the superscripts denote specific iterations. Suppose that at some iteration at least one constraint is 
violated by the point Let this constraint be the jth hyperplane. 


ai Xj -t- ^2 a 2 i- ... -r a^j^ x^^ d. . 


(28) 


Substituting (27) into (28) gives 


ajjCxii + aiHivf(xi)i) + a2Hx2^ + aitf Vf(x»)2) + - + aj(x^i + HiVf(xi)n,) = bj , 


where Vf(x^)j. is the kth element of the search vector. Upon gathering terms 


(a) )Txi + Vf (x^) = bj . 


Then 

Vf(xi) = bj - (ai)Txi . 

Since (aj)"*"rfVf(x^) is a scalar, can be computed as 
ai = (bj - (ai)Txi)/((aj)%i Vf(xi)) . 
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This is used to characterize the distance from to a point at the intersection of the search vector and 
each constraint hyperplane that is violated by In the case where x^ is an interior point in feasible space 
and is an exterior point, the smallest would be selected. The selected will be denoted as a*. The 
knowledge of will be used to eliminate the nonbinding constraints along a search vector. Figure 6 shows 
graphically that the search vector p^ from x^ may intersect many constraints. Note that g^ should be the 
only binding constraint and also has the smallest 



Figure 6. Intersection of search vector and constraints. 

The method of calculating the projection matrix was developed from the well-known Gram-Schmidt 
technique [15]. All the Vn^(x)’s of the binding constraints span the q-dimensional space Q. From abstract 

algebra it is known that the basis of Q is not unique and by employing the Gram-Schmidt method an ortho- 
normal basis for Q can be derived. Choosing ^ as an arbitrary element of the constraint manifold, other 

elements can be found by the method used by A. O. Morris [16]; 


i^ = Vnj(x) 

i^ = Vnj.(.i(x)- (Vnj+i (x)’^ i^) iVHi^ 

P = Vnj+2(x) - (Vnj+2(x)T i^) iV Hi* 

i^l = Vnq(x) - (Vnq(x)^ i*) i*/lli* iP - ... - (Vnq(x)^ i^l'^/lli^l"* 11^ . 


I|2 

l|2-(Vnj+2(x)'^i^) i^/lli^l|2 
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The P vectors are changed to a unit vector by replacing each P by its unit vector where 


uJ = ij/llij|| . 


The projection of V f(x), where x is a point in the intersection of the q constraint hyperplanes, onto this 
vector set will be 


p = (I - u^(u^)^ - u^(u^)"^ - ... - u*^(u*^)"^) V f(x) . 


Let 


q 

A = ^ < 

j=l 

Then 


p = (I-A)Vf(x) . 


The method of calculating this A matrix is crucial to the efficiency of the method. From (29) it can be seen 
that the Gram-Schmidt method is also a projection of each succeeding Vnj(x) upon the previous vector set. 

Using this, the following recursive method was developed to update A: 


ii+1 =(I- A|)Vni+l(x) 


(30) 


This vector is normalized to a unit vectors as before. Then 


Aj.^j = Aj + (u^"*"^)"^ • 


(31) 


When updating A by the addition of constraints, they can be added to the A matrix in any order. The 
removal of constraints is not so simple. In the algorithm, all elements of the A matrix are set to zero and 
are then recalculated. Some of the previous computational work is saved which allows a portion of the A 
matrix to be recalculated with very little computation. Each iteration that updates the projection matrix has 
the same properties as the Gram-Schmidt method in that it is sensitive to the order in which the vectors are 
used in the algorithm. Because of this, the algorithm does not produce a unique projection matrix at inter- 
mediate steps. To eliminate a constraint, simple subtraction cannot be used. The vectors must be removed 
in the reverse order in which they are added. This continues until all the nonbinding constraints have been 
removed. This process will also remove some binding constraints which must then be re-added. The same 
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number of operations are required to remove a vector as to add a vector. If the vector to be removed had 
been added at the mid-iteration of the algorithm then, as many operations would be required to remove the 
constraints in reverse order and then regenerate the matrix, as generating the entire projection matrix would 
require. To minimize the computational effort to regenerate the projection matrix, the program presented 
saves the u’s of the appropriate binding constraints. By zeroing all elements in the projection matrix and 
updating it with only the active constraints, the effect is the same as that of dropping aU the nonbinding 
constraints at once and adding only the binding constraints. 


VI. CONCLUSIONS 


The theoretical development for the method of nonlinear optimization presented in Section IV 
shows the applicabihty of this method and some of the advantages over the Rosen method. The method 
proposed by this paper provides similar results with fewer computations per iteration. It has been noted 

T 1 

earlier that the method presented here eliminates the need to calculate (N^N)' as required by the Rosen 
method. This not only saves computer time, but avoids a potentially formidable matrix inversion problem 
as the system becomes large. Another advantage of the method proposed here is that there is no limit to the 
number of constraints that can be simultaneously added or removed from the binding set. The Rosen 
method is limited to one addition or removal at a time. 

A computer program was developed which implements this method. It is constructed in a modular 
format to facilitate understanding and provide efficiency in application. The program was applied to a 
variety of test problems and was successful in finding the optimum in each case. These cases were chosen to 
evaluate the program’s ability to converge while encountering some of the difficulties that are inherent in 
projection methods. A good comparison between this method and the one used by Rosen could not be made 
because of possible differences in the type of computer used and the programming structure. The program’s 
computer run time is strongly dependent on the techniques used in the line search and the speed and register 
capacity of the computer. The cases were chosen not as a comparison to problems worked by other pro- 
grams, but to demonstrate the abilities of this program. The only direct comparison is based on the number 
of computations required by each method to update the projection matrix. This should provide some idea 
of the relative speed of operation during the calculation of the projection matrix. Table 1 in Section I shows 
that the method proposed by this paper should be more efficient with computer time when calculating the 
projection matrix. 

An area of future research would be an investigation of the differences of the projection matrix used 
by the two methods. A preliminary look at the projection matrix generated by each method for several 
cases showed that they were the same. If the two methods do generate the same projection matrix for all 
cases, this should be proven rigorously. 
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APPENDIX A. USERS GUIDE 


The computer program for the solution method developed in this paper requires two subroutines 
and a data deck to be supplied by the user. The two subroutines are incorporated into the program as sub- 
routines called FUNCT and GRAD. They take the following form: 


SUBROUTINE FUNCT(X,F,N1) 
DIMENSION X(N1) 

F = Objective function 

RETURN 

END 

SUBROUTINE GRAD(X,G,N1 ,NUM) 
DIMENSION X(N 1 ),G(N 1 ) 

G(l)=the first gradient element (3f(x)/c»xj) 


G(Nl)=the Nlth gradient element (3f(x)/3xj^j) 

DO 1 1=1, N1 
1 G(I) = NUM*G(I) 

RETURN 

END 

where 

N1 is the number of independent variables 

F is the function to be extremized 

X(i) is the ith variable element 

G(i) is the ith gradient element of F 

NUM=1 for maximization and = -1 for minimization. 
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The coefficients and right-hand side constrant value for the system of linear constraints is read into 
the program from the data deck with the initializing data. The description of each of the data cards follows; 

DATA CARD 1 : 

(MAX/MIN),N1 ,NN,NC,EP1 ,EP2,NLEC,NEC,NGEC,NGEX 
These are read according to FORMAT(A3,3I3,2E 14.7,413) where 
MAX for maximization or MIN for minimization problems 


N1 — Number of variables in the problem 

NN — Maximum number of unconstrained iterations allowed. This wiU limit the number of 

iterations in case of slow convergence to the solution for the unconstrained DFP method. 


NC - NLEC+NEC+NGEC+NGEX 

EPl — If lF(x^‘*'^)-F(x^)l < EPl; EPl is typically = 0.0001 and detects small changes in the 
function value. This is a condition for convergence. 

EP2 — If II VF(x^‘*’^)|| < EP2; EP2 is typically equal to 0.0001 and is the norm of the function 
gradient. It is also used as a condition for convergence. 

NLEC — The number of < constraints 

NEC — The number of = constraints 

NGEC — The number of > constraints (excluding NGEX) 

NGEX — The number of x^ > 0 i=l,2,...,n constraints. 


DATA CARD 2 : 

X0(I),I=1,N1 

XO is the initial point and is read according to FORMAT(10E13.6). If NC=0 (unconstrained problem), the 
following data cards are omitted. 


DATA CARD 3 : 

BC(I),I=1,NC 

BC(I) provides the right-hand side of the constraints. The BC(I) values must be provided in the following 
order (1) less than or equal to constraints, (2) equality constraints, (3) greater than or equal to constraints, 
and (4) nonnegative variable constraints. This ordering is illustrated by the following example: 
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gi(x)<b i= 1,2,...,NLEC 
gj(x) = b j= 1,2,.., NEC 
gk(x)>b k= 1,2,...,NGEC 
gg(x) = xg > 0 

The BCs represent the b values and are read according to FORMAT(10E13.6). The problem requires a 
standard formulation such that x^>0 for all i. 

DATA CARDS FOR CONSTRAINTS : 

The coefficients for each of the constraints are read according to FORMAT(10E13.6). The order for 
the constraints must match the order for the right-hand side constants on data card 3. All coefficients must 
be included, i.e., N1 coefficients for each constraint including zeroes. 

Example problem and corresponding user input: 

Consider the following example : 

max F = (xj - 3)2 + 9)X2 - 5)2 

X0 = (0.5,-1.2,) 

Such that 

2xj + X2 < 20 

xi + 2x2 < 40 
xj + 2x2 < 30 
9 xj + 6x2 =100 
xj > 0 
X2>0 . 

SUBROUTINES : 

SUBROUTINE FUNCT(X,F,N1) 

DIMENSION X(N1) 

F=(X(I)-3.)**2+9.*(x(2)-5.)**2 

RETURN 

END 
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SUBROUTINE GRAD(X,G,N 1 ,NUM) 
DIMENSION X(N 1 ),G(N 1 ) 
G(l)=2*(X(l)-3.) 
G(2)=18.*(X(2)-5.) 

DO 1 1=1, N1 
1 G(I)=NUM*G(I) 

RETURN 

END 

DATA CARDS 

MAX 2 70 6 0.0001 

0. 5 -1.2 

20. 40. 30. 100. 

2 . 1 . 

1 . 2 . 

1 . 2 . 

9. 6. 

1 . 0 . 

0 . 1 . 


0.0001 3 1 0 2 
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APPENDIX B. TEST PROBLEMS 


This section presents the results of several test problems that were solved with the technique 
developed in this paper. Unconstrained problems were selected to test the DFP method. Linearly constrained 
problems were then selected to test the projection technique. The program was initialized at infeasible points 
for some cases to test the algorithm that generates a feasible starting point. 


UNCONSTRAINED TEST PROBLEMS AND RESULTS 

The following optimization problems were unconstrained. The program uses the DFP method with 
the approximate Hessian inverse matrix reset to I after N+1 steps. The restrictions for x^ > 0 are relaxes for 
the unconstrained problems in which NC = 0. 

T1 : BANANA FUNCTION 


N-1 

MIN f(x) = Y, [ 100(Xi+i - Xi2)2 + (1 _ x-)2] 
i=I 


The starting point with N = 3 is x° = (-l.,0.,-l.). The minimum was found after 73 iterations of the DFP 
method. The program execution time was O.1 125 min on a SIGMA-V computer. The minimum point of the 
objective occurred at x* = (l.,l.,l.) with f(x*) = 0.0. This problem is also known as Rosenbrock’s Function 
when N = 2. The Banana function is considered to be a severe test for nonlinear optimization algorithms due 
to the steep valley generated by the coefficient value of 1 00. 


N-1 

T2: MIN f(x) = Y n0(xj+j - Xj2)2 + (i _ x^2] 
i=l 


The starting point with N = 3 is x° = (6.67,6.67,0.). The minimum of the objective occurred at x* = 
(l.,l.,I.) with f(x*) = 0.0. This problem is similar to Tl. The difference is the reduced steepness of the 

valleys formed by the terms. Computer time for execution was 0.0284 minutes and the number 

of iterations was reduced to only 12. 

T3: MIN f(x) = (xj - 3)^ + 9(x2 - 5)^ 

The starting point is x° = (l.,I.) and the minimum occurred at x* = (3. ,5.) with f(x*) = 0.0. Three iterations 
of the DFP method were performed. The computer time for execution was 0.0159 minutes. 

10 

T4: MIN f(x) = X) (x - 1 0)^ 
i=l 
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The starting point is x° = (0.,1.,2.,3.,4.,5.,6.,4.,2.,4.). The minimum occurred at x* = (10., 10., 10., 10., 10., 
10.,10.,10.,10.,10.) with f(x*) = 0.0. Only one iteration of the DFP method was required since this is a 
spherical function. The total execution time was 0,0292 minutes. 


CONSTRAINED TEST PROBLEMS AND RESULTS 

For these problems the iterations when using a projected search vector are separated from those of 
the DFP method. The Gradient Projection Method makes no use of the Hessian (or its inverse). For the 
constrained problems, once the solution sequence encountered the boundary of feasible space, it remained 
on the constraint manifold. All these examples were initialized outside of feasible space in order to test the 
routine that finds a feasible initial point. 


2 

T5 : MIN f(x) = ^ [ 100(xj+i - Xj2)2 + (i _ x^2] 
i=l 


subject to 

2xj + X 2 + x^ < 20 

X 2 + 2x2 + 4 x 3 < 40 
xj + 2x2 + 2 x 3 < 30 
9 xj + 6x3 X3 == 100 
lOxj + 20x2 + X3 > 100 

xj > 0 , i = 1 ,2,3 


The minimum of the objective occurred at x* = (6.67,6.67,0.) with f(x*) = .33E+6. The search terminated 
at this point due to satisfying a program convergence test. Since it is not possible to drop equality con- 
straints from the projection matrix, the last point is taken as a constrained extremum. The minimum was 
found with two projection iterations and no DFP iterations. The total execution time was 0.0630 minutes. 

2 

T6: MAX f(x) = lOOfx^+j - + (1 - x^^] 

i=l 


The constraints are the same as those in T5. The maximum was found at x* = (1.67,14.17,0.) with f(x*) = 
.34E+6. At this point a constrained maximum occurred caused by the restriction of not dropping the 
equality constraints from the projection matrix. This demonstrates the abiUty to move around on the 
equality constraint manifold in the search for an extremum. This point was found after three projection 
iterations and no DFP searches. The execution time was 0.1 186 minutes. 
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T7: Max[Tl] 


such that 

2xj + X2 + 4xg < 20 
X| + 2x2 ^ 

xi + 2x2 + 2 x 3 <30 

9xj + X2 + X3 < 100 
lOxj + 20x2 + X3 > 100 
X| > 0 , i = 1,2,3 


The constrained maximum was found at x* = (0.,15.,0.) with the value of f(x*) = .65E+5. The maximum 
was obtained in one projection search, and the program execution time was 0.0703 minutes. 

The results for T5, T6, and T7 are shown in Figure B-1. The tests T5 and T6 have an equality con- 
straint that the extremum is required to satisfy. The equality constraint is constraint number 4. The test T7 
had no equality constraints and the extrema for this problem was found at a vertex. The equality constraint 
was changed to an inequality. The test T5 is a minimization and the test T6 is a maximization of the objec- 
tive function which explains why the extrema are found as far away from each other as possible and still 
lie on the equality constraint. 

T8: MIN[T7] 

The constrained minimum was found at x* = (2. 8, 3. 5, 2. 7) with f(x*) = .lE+5. The minimum was found 
after eleven projection iterations. The execution time was 0.1020 minutes. Constraints 1 and 5 are binding 

as shown in Figure B-2. The initial point was x^ = (0.,5.,0.). From this point the search proceeded along 
constraint 5 until constraint 1 was encountered. The minimum of the objective function was found in the 
seam formed by constraints 1 and 5. 

T9 : MIN [f(x) = 5e^^^^ + X2X4^ - X3 sin(2(x2)) + x^^ - lOxgX^j 
such that 

X 2 + 10 x 4 + xg < 35 

X3-5X5 + Xg-X 7.<8 

xi + X 3 + Xg < 10 

X j + 3x^ + Xg + 7xy < 200 

-X 2 + 1 5 x 3 + Xg - X 7 + 3xg = 25 
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X4 - Xs + 5xg - lOxg = 27 
-3x2 ~ ^^3 ^5 ^ 

X]^ + 15x2 + X 5 + 8 xg >17 

xj - 2x2 - 8xg + 5x4 + X5 - 7xg + lOxj + 3xg > 50 

X|>0 , 1=1,2,..., 8 

The starting point x° = (-1 .,0.,-l .,0.,-l .,0.,-l .,0.) was outside the feasible region and the program generated 

a new starting point x° = (0.,0., 1.89014,3.5 ,0.,4.7,8.0521 1,0.), The minimum of the objective function 
occurred at x* = (0.,0., 1.89014,3 .5,0. ,4.7, 8.0521 1,0.) with f(x*) = 5.727. This constrained minimum was 
found after six projection iterations. The extremum was found in constraints 1, 5, 6, 9, 10, 11, 14, and 17. 
The execution time was 0.1224 minutes. This is a vertex of the constraints. 

TIP : MIN [f(x) = (xj - 18)2 + (X2~ 18)2] 

such that 

-20xj + 4x2 < 20 

-4.7x2 + 2x2 < 14 

-1.5x2 + X 2 < 1 1 

-2x2 2^2 ^ 20 

0.4x2 + 2x2 < 40 

2x2 + 0.4x2 ^ 42 

x->0 , i = 1 , 2. 

The starting point was x° = (0,0) and the minimum point was x* = (17.54,16.49) which is a vertex. Figure 
B-3 shows the search path the program followed to acquire the constrained extremum. This shows the 
advantage of the nonlinear method over the Simplex method for some problems in that it can traverse 
feasible space. A Simplex method would have taken five more iterations. The execution time was 0.0109 
minutes. 
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APPENDIX C. PROGRAM DESCRIPTION 


Appendix C contains a brief description of the function of each subroutine of the program and how 
they interact. The description of each subroutine generally follows the order that data flows through them. 
In some cases the logic flow in the program is controlled by flags that are passed between subroutines. The 
program is modular with each of the modules or subroutines generally performing one logical function. 
Examples of this are found in the MAIN program where the DFP method is located or in subroutine SEEK 
where the line searches are performed. 


MAIN PROGRAM 

The main program initializes the parameters, which include the number and type of each constraint, 
and assesses the supplied initial starting point to determine if it is in feasible space. If it is not or if the 
problem includes any equality constraints, the MAIN program calls SUBROUTINE FEASPT to internally 
generate an initial starting point in feasible space. The program then calls SUBROUTINE FUNCT and 
SUBROUTINE GRAD to get initial values for f(x) and V f(x). The Hessian is initialized to I and SUBROU- 
TINE SEEK is called. For the case where the extremum of the function is located in feasible space or where 

there are no constraints, then SEEK returns new values for x^, f(x^), and V f(x^) that represent an approxi- 
mation to the line search extremum. These values are used to calculate an update Hessian matrix. The pro- 
gram tests the new point against the criterion for a local extremum. If it has not satisfied this criterion 
another sesu’ch direction is generated and another line search is made. This continues until a local extremum 
is found or until the maximum iteration Umit is acquired. In the case where the extremum is not located in 
feasible space, the constrained extremum point will be returned to MAIN from SUBROUTINE CHK for 
printing and program exit. 


SUBROUTINE FUNCT and SUBROUTINE GRAD 

These are discussed together since they are both user supplied subroutines. The equation for f(x) 
is supphed by the user in FUNCT, and the gradient functions Vf(x) are given in GRAD. The gradient vector 
returned from GRAD has the proper sign for maximization or minimization. 


SUBROUTINE SEEK 

This subroutine calculates the DFP search vector when the search is being conducted from the 
current search point which is unconstrained. When constraints are determined to have been violated, the 
search vector is supplied by SUBROUTINE CHK. The program is sensitive to the magnitude of the search 
vector. It adjusts the initial step sizes along this vector inversely to the magnitude of the vector. As it steps 
along this vector, new values of f(x) are calculated. When a change of direction from decreasing to increasing 
values of f(x) along the search direction is detected (extremum is bracketed), a quadratic curve fit of the 
last three points is made in order to estimate the location of the extremum. Figure C-1 shows the curve fit 

made through y(a) = f(x^ + oi^H^Vf) point for successive Oi’s. This method estimates a* for the, minimum 

point of a quadratic curve fit where f(x^'*‘^) = min[f(x^ + ce^H^Vf(x^)]. This approximation of the actual 
extremum gets better near convergence. In the case where the change of direction in the values of f(x) 
occurs on the first step along the search vector, the quadratic curve fit is made based on the original starting 
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y(a) 



Figure C-1. Extremum estimation along search vector. 

point of the search vector, the gradient at that point, and the first step point. SUBROUTINE SEEK takes 
ten steps along the search vector. If the values of f(x) have not changed direction, then an error flag is 
returned to the MAIN PROGRAM. If the problem is of such a complexity that a reasonable estimate cannot 
be made, then SUBROUTINE FIB is called to supply this point. SUBROUTINE CHK is always called to 
determine if the point found is in feasible space. If it is not, then CHK returns a new search vector. A return 
is made to MAIN if a constrained extremum was determined. 


SUBROUTINE FIB 

This subroutine performs an eight-step Fibonacci Search along the vector supplied by SEEK. The 
new point which is returned is an approximation of the extremum point along the search vector. An eight- 
step Fibonacci search should be sufficient to provide the proper point within a small error. This type of 
search method is used since the number of steps have been specified in order to obtain the minimum interval 
of uncertainty for the placement of eight points. 


SUBROUTINE CHK 

This subroutine tests for constraint violations. If no violations have occurred, it returns the same 
search vector to SEEK. When new constraints are violated by a search point, CHK selects those that are 
violated. These are added to the projection matrix. If the number of these constraints equals or exceeds the 
number of dimensions, then SUBROUTINE DROP is called to determine which constraints are to be 
dropped. SUBROUTINE CHK determines when a constrained extremum has been encountered. When either 
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(1) the projected search vector magnitude becomes smaller than a user set tolerance or (2) the search vector 
points away from any equality constraint, a constrained extremum is assumed to have been found. When 
SUBROUTINE DROP returns a new search vector and a new point, SUBROUTINE CHK determines when 
all new constraints have been acquired and then searches these to see if all of the equality constraints are 
contained in the projection matrix. If they are not, then the original point sent to DROP is assumed to be 
the constrained extremum. If the same constraints are acquired again, then this point is assumed to be con- 
strained extremum. The projection matrix is returned from SUBROUTINE SPAN. It uses the set of con- 
straints supplied by CHK. 

SUBROUTINE SPAN 

This subroutine uses an algorithm that is based on the Gram-Schmidt method to calculate the pro- 
jection matrix. Any number of new constraints may be added without reinitializing the matrix. When con- 
straints are dropped, all elements in the projection matrix are set to zero and it is recalculated. 


SUBROUTINE DROP 

In this subroutine the search vector is tested to determine if it points into or out of feasible space. 
If it points into feasible space, then only the equality constraints are kept to calculate the new projection 
matrix. If there are no equality constraints, then a flag is set which will cause the program to return to 
MAIN and initiate the DFP Search method from this point. If the search vector points out of feasible space, 
then DROP calculates a new starting point at a very small distance from the constraints and sets the search 
vector equal to the gradient direction. A return is made to CHK to continue the search. 


SUBROUTINE FEASPT 

This subroutine generates an initial feasible point using Phase I of the two-phase Simplex method. 
The Phase I method finds a point that lies on the boundary of feasible space which satisfies all of the 
equality constraints. This point is transferred to the MAIN Program. If Subroutine FEASPT fails to find a 
feasible point or if more equality constraints are present than the number of variables, a message is printed 
indicating the type of difficulty encountered, and an error flag is set that terminates the program. When too 
many equality constraints are present, it indicates that the system has linear dependency, redundancy, or 
is an inconsistent set. 
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APPENDIX D. PROGRAM LISTING 



OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO 


4 ^ 

o\ 


it* X- *************** *****■****•********)<■*■)(•■)(•****■)*•*************** 

THE CONV'jMl’ION WILL 3E 'Mk'f THE NORMALS TO TdE 
COx'JSTRAINT PLANES WILL 3E NEGATIVE EOR THOSE THAT 
POINT INTO THE .FEASIBLE SPACE AND POSITIVE .FOR THOSE 
THAT POINT INTO IxNEEASIBLS SPACE. 

X-*** ******************** ****it **************** **********x**** 

THE CONSTRAINTS WILL BE ORDERED IN THE PROGRAM IN 
THE FOLLOWING WAY; 

NO. OF AX<B CONSTRAINTS = NL.3C 

NO. OF AX=;B constraints = NEC 

NO. OF AX>B CONSTRAINTS = NGEC 

THE X>0 CONSTRAINTS WILL NOT BE IN NGEC 

NO. OF X>0 CONSTRAINTS = NGEX 

************ X**************** X*** ********* *****X X******* X*** 


N1 

NN 

EP1 

EP2 

Y = 

YL 

GR1 


NO. OF DIMENSIONS 
MAX. NO. OF ITERATIONS 

>1U(N+1 )-X(N)11THIS is for convergence CRITERIA. 
>ltGRAD(N)![ THIS IS A COxIVERGE.NCE CRITERIA. EP2=. 
OBJ. F'JNCT. VALUE 
LAST OBJ. FUNCT. VALUE 


EP1 = .0001 GEN.ERALLY 
0001 GE.NERALLY 


= GRAD. 


GRO = 
FLAG1 

FLAG2 

FLAG3 

FLAG4 

FLAG5 


LAST GRAD. 

= 1 SPECIFIES THAT THE) 
AND IS COMPLETE. IT 
THAT Txlfi 


= 0 


'H' MATRIX HAS BEEN UPDATED N1 
WILL BE USED ONE MORE TIME. 
'H' MATRIX IS NOT YET COMPLETE 


TIMES 


A CONSTRAINT HAS BEEN 
IN PLANE VIOLATED. 


1 USED IN SEEK TO ALERT IT THAT 
VIOLATED AND DOES A LINE SEARCH 

2 CAUSES SEEK TO SET FLAG4 = 1 

1 SPECIFIES THAT TflPJ COxIDITIONS FOR A CONSTRAINED 
MAX (MIN) HAS BEEN MET. 

2 THE I'lAX. NO. OF TRIES AT PINDIlNC 


AN EXTREMA ALONG 


THE SEARCH VECTOR HAS BEEN DONE. 


= ALERTS PROGRAM THAT IT WAS IN A CONSTRAINT BUT IS 
NOW IN FEASIBLE SPACE. RESETS H=I,B.%C=0 AND AAT2=0 
= ALERTS PROGRAM THAT TH.fiRE IS LINEAR DEPENDENCE AMONG 



ooooooooo 


THE EQHALITY CONSTRAINTS. 

KK = THIS IS A COUNT OP THE TOTAL NO. OP ITERATIONS THROUGH THE 
HESSIAN. IT IS USED TO PORCE THE PROG. THROUGH AT LEAST 
N1 UPDATES OP 'H'. USED TO ABORT PROG. IP ERROR BOUNDS 
CAN'T BE MET BY SPECIFIED NO. OP ITERATIONS. 

E = RESETS H = I APTER N1 ITERATIONS. 

NC = NO. OP CONSTRAINTS INCLUDING THE Xi>0 CONSTRAINTS 
NUM = THE VARIABLE THAT DESIGNATES EITHER MAX OR MIN. 

IT DOES THIS BY IT'S SIGN (+/ IS MAX, -/ IS MIN) 

COMMON XOdOKXI (10),H(10,10),A(20,10),GR0(10),BC(20), 

1 N 1 , NC , AAT2 ( 1 0 , 1 0 ) , KOP , NC VO ( 20 ) , NLEC , NEC , NGEC , NGEX , KC 
DIMENSION GR1 00),XD(10),GRD(10),P(10), 

1B(10,10),B1 (l0),C(10,10),Cri(l0),CHN(lO,10) 

INTEGER PLAG1 , PLAG2 , PLAG3 , PLAG4 , PLAG5 , NC , TYPE , QMIN , QMAX 
DATA QMAX, QMIN /4HMAX ,4HMIN / 

Q***)t* K0P=0 MEANS THAT THE SEARCH MUST START IN FEASIBLE SPACE. 

C INITIALIZE 
PLAG5=0 
K0P=0 
KK=0 
XC=0 
K=0 
NUM=1 

READ(5,300)TYPE,N1 ,NN,NC,EP1 ,EP2, NLEC, NEC, NGEC, NGEX 
WRITE(6,320)TYPE,N1 ,ilN,NC,EP1 ,£P2, NLEC, NEC, NGEC 

320 P0RMAT('NE ARE LOOKING FOR A ',A3/ 

1 'NO. OP DIMENSIONS =’,I3/’MAX. NO. OP ITERATIONS =',13/ 

2 'NO. OP COilSTRAINTS =’,I3/'EP1 & EP2 =',2E13.6/ 

3'NO. OP LE CONSTRAINTS=' ,I3/'N0.0P = C0NSTRAINT3= ' , 13/ 

4'NO. OP GE CONSTRAINTS=' ,13) 

READ(5,301 )(X0(I),I=1 ,N1 ) 

WRITE(6,321 )N1 , (X0(l),I=1 ,N1 ) 

321 P0RMAT('X0(I),I=1 ,N1 '/nE13.6) 

IP(NC.EQ.O)GO TO 33 

DO 34 1=1 ,NC 
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4!>. 

00 


34 READ(5, 301 )BC(T ) 

DO 1 T=1,NC 

1 READ(5, 301)A(T, J) ,J=1 ,M1 
WRTTE(6,322) 

322 FORMATCTHTS IS THE "A" MATRIX') 

NO = NC-NGEX 

WRT.TE(6, 318)(N1 , ( A{I, J) ,J = 1 ,N1 ) ,1 = 1 ,N0) 

318 F0RMAT(NE1 3.6) 

WRITE(6,400)NO, (BC(I ),I = 1 ,N0) 

400 FORMATCBCr ',WE13.6) 

***** IF THERE ARE EQUALITY CONSTRAINTS THEN THE PHASE 1 METHOD 
WILL BE USED AS THIS WILL PUT THE INITIAL POINT ON THE 
EQUALITY CONSTRAINTS AS REQUIRED. 

IF(NEC.NE.O)GO TO 32 
DO 29 I=1,NC 
BC1=0.0 
DO 28 J = 1 ,N1 

28 BC1=A(I,J)*X0(J)+BC1 
IFd.GT.NLEOGO TO 40 
IF((BC(I)-BC1).LT.O.O)GO TO 32 
GO TO 29 

40 IF(I.GT.(NLEC+NEC))GO TO 41 
IF(ABS(BC(I)-BC1).GT. 1.E-12)G0 TO 32 
GO TO 29 

41 IF((BC(I)~BC1).GT.O.O)GO TO 32 

29 CONTINUE 
GO TO 33 

32 CONTINUE 

***** THIS ROUTINE IS THE SIMPLEX PHASE 1 METHOD. IT PIVOTS 

THROUGH THE VERTICIES UNTIL IT LOCATES ONE IN FEASIBLE SPACE. 
CALL FEASPT(FLAG5) 

IF(FLAG5.EQ.1 )G0 TO 36 

C***** IT IS ASSUMED THAT THE CONSTRAINTS ARE LINEARLY INDEPENDENT. 
IF(N1.GT.NEC)G0 TO 33 
DO 35 1=1 ,N1 



o o o 


35 XUT) = XO(.T) 

GO TO 23 
33 CONTINUE 

IF(QMAX.EQ.TYPE)NUMr1 
IF(QMIN.EQ.TYPE)NUM=-1 
323 FORMAT('TH£ DATA HAS BEEN READ') 

C**it**xO = INITIAL POINT 

C*i.**»Y = VALUE OF THE OBJ. FUNCT. 

C*»***N1 = NO. OF DIMENSIONS 
= NO. OF CONSTRAINTS 
CALL FUMCT(XO,Y,N1 ) 

210 CONTINUE 
FLAG3=0 
YL=Y 

350 F0RMAT('X0(I),Y'/NEU.7) 

C GRAD. OF THE FUNCT. AT THE POINT XO 
CALL GRAD(X0,GR0,N1 ,NUM) 

351 FORMAT( 'GRO(I)' /NE14.7) 

100 CONTINUE 

C*****THIS SETS THE HESSIAN = I 
DO 3 1=1, N1 
DO 2 J=1,N1 
AAT2(I,J)=0.0 
H(I, J)=0.0 

2 IF(I.EQ.J)H(I,J)=1.0 
NCVO(I)=0 

3 CONTINUE 
FLAGUO 
FLA34=0 
K=0 
YL=Y 

309 CONTINUE 

***»*THIS IS THE LINE SEARCH MAX(MIN )f(X1+^HGRAD) 


- 1 ^ 

VO 
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CALL SEEK (Y,YL, FLAG1 , FLAGS, FLAG4,NUM,RH0,P) 

C 

C***** IF FLAG3=2 THE MAX NO. OF ATTEMPTS HAS BEEN DONE ALONG 
the search VECTOR. 

C 

TF(FLAG3.EQ.2)GO TO 27 

IF FLAG3=1 A CONSTRAINED MAX(MIN) HAS BEEN FOUND 
IF(FLAG3.EQ. 1 )G0 TO 23 

***** IF FLAGi| = 1 THE GRAD NOW POINTS BACK INTO FEASIBLE SPACE 
IF(FLAGi4.EQ. 1 )GO TO 200 

***** IF FLAG1=1 A LINE SEARCH WITH THE COMPLETE HESS HAS BEEN 
COMPLETED AND NOW RESET HESS=T 

IF(FLAG1.EQ. 1 )GO TO 22 

CALL GRAD(X1,GR1,N1,NUM) 

WRITE(6,3U) 

314 F0RMAT(5X, ’XI ' , 1 2X , • XO ' , 1 2X , ' GR 1 • , 1 2X , ' GRO ' ) 

WRITE(6,315) (X1(T) ,X0(T),GR1(I) ,GR0(I) ),I=1 ,N1 

315 F0RMAT(4E14.7) 

DO 4 1=1, N1 
XD(I)=X1(I)-X0(I) 

4 GRD(I)=(GR1 (I)-GR0(T))*NUM 
B2=0.0 
C2=0.0 
DO 20 1=1 ,N1 
B2=B2+XD(I) **2 
20 C2=C2+GRD(I )**2 
B2=SQRT(B2) 

C2=SQRT(C2) 

GO TO 22 



200 CONTINUE 

IF(FLAG1.EQ. 1 )G0 TO 19 

c*»***resets the b&c matrices = 0 

K0F=0 

DO 16 I=1,N1 
NCV0(I)=0 
CH(I)=0.0 
DO 15 J=1,N1 
AAT2(I, J)=0.0 
B(I,J)=0.0 

15 C(I,J)=0.0 

16 CONTINUE 

IF(FLAG4.EQ. 1 )G0 TO 100 

B3=0.0 

DO 5 1=1, N1 

5 B3=B3+P(I)*GRD(I ) 
B3=ABS(RH0)/B3 

C 306 F0RMAT('B3=' ,E14.7) 

DO 6 I=1,N1 
DO 6 J=I,N1 
B(I,J)=P(I)*P(J)*B3 

6 B(J,I)=B(I,J) 

DO 7 1=1, N1 
CH(I)rO. 

DO 7 J=1,N1 

7 CH(I)=CH(I)+H(I, J)*GRD(J) 
CHD=0. 

DO 8 I=1,N1 

8 CHD=CHD+GRD(I)*CH(I ) 

DO 9 1=1, N1 

DO 9 J=1,N1 

9 CHN(I, J)=CH(I)»GRD(J) 

DO 11 1=1, N1 

DO 11 J=1,N1 
C(I,J)=0. 
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K> 


DO 10 J1=1 ,N1 

10 C(I, J)=C(I,J)+CHN(I, J1 )*ii(J1 ,J) 

11 C(I,J)=G(I, J)/CHD 
DO 12 1 = 1 ,i11 

DO 12 J = 1 ,N1 

12 d(I, J)=ri(I, J)+B(I,J)-C(I,J) 

19 CONTINUE 

KK=KK+1 

K=K+1 

DO 21 1 = 1 ,N1 
X0(I)=X1 (I) 

21 (JR0(I)=GR1 (I) 

XL=Y 

GO TO 309 

250 CONTINUE 

DO 251 1=1, N1 
X0(I)=X1 (I) 

251 CONTINUE 
GO TO 210 

22 CONTINUE 

■****MUST GO THROUGH ^^T LEAST N1 CALC'S OE 'H' FIRST TIME 
IP(KK.LT.N1 )G0 TO 24 


♦****THIS CHECK IS TO DETERMINE IP YOU'RE CLOSE ENOUGxH 


IP(B2.LT.EP1 .AND.C2.lt. EP2)G0 TO 23 
24 CONTINUE 

IF ( FLAG1 . EQ . 1 ) GO TO 1 00 

IF(K.EQ.N1 )PLAG1=1 

*****IP KK=HAX NO. OP ITERATIONS THEN STOP 


IP(KK.EQ.NN)G0 TO 23 



c 

C*k»*«tHIS allows one step with the completed 'H' matrix 

c 

GO TO 200 
23 CONTINUE 

IF(FLAG3.E0.0)G0 TO 27 
WRITE(6,305)TYPE 

305 FORMATCA CONSTRAINED ',A3,' HAS BEEN FOUND.*) 

27 CONTINUE 

IF(FLAG3.EQ.0)G0 TO 37 
IF(FLAG3.NE.2)GO TO 13 
WRITE(6,32^) 

324 FORMATCTHE line SEARCH HAS TAKEN THE MAX. NO. OF STEPS') 
CALL GRAD(X1,GR1,N1,NUM) 

GO TO 37 
13 CONTINUE 
DO 38 I=1,N1 
X1(I)=X0(I) 

38 GR1(I)rGR0(I) 

37 CONTINUE 

WRITE(6,303)KK,KC 

303 FORMATCNO. of unconstrained ITERATIONSr* ,14// 

*'N0. OF CONSTRAINED ITERATIONS: 14 ) 

WRITE(6,304)Y,N1,(X1(I),I = 1 ,N1 ),N1,(GR1(I),I = 1 ,M1 ) 

304 FORMAT( 'F(X*)=' ,E1 4.7//'X*=* ,NE1 3.6// 

‘•'GRADIENT F(X *) = ' , NE 1 3. 6 ) 

300 F0RMAT(A3,3T.3,2E14.7,4I3) 

301 F0RMAT(10E13. 6) 

36 CONTINUE 

END 



SUBROUTINE FUNCT(A,B,N) 

DIMENSION A(N) 

B=100.*(A(2)-A(1 )**2)**2+(1.-A(1 ))**2+ 
*100.*(A(3)-A(2)*«'2)**2+( 1.-A(2))»*2 
RETURN 
END 


SUBROUTINE GRAD(A, B, N, NUM) 

DIMENSION A(N),B(N) 

B(1 ) = -40.»(A(2)-A(1)**2)*A(1 )-2.*(1 .-A(1 )) 

B(2)=20.*(A(2)-A(1 )*'*2)~2.*(1.-A(2) )-4 0.*( A(3)-A(2)*^*2)*A(2) 
B(3)=20.*(A(3)-A(2)**2) 

IF(NUM.GT.O )RETURN 
DO 1 I=1,N 
1 B(I )rNUM*B(I ) 

RETURN 

END 
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SUBROUTINE SEEK(Y , YL , FLAG1 , FLAG3, FLAQil , NUM , RHOM , DD ) 

B = B MATRIX OF CONSTRATNTS(A ' X=B) 

NC = NO. OF CONSTRAINTS 

A = MATRIX OF THE CONSTRAINTS COEFFICIENTS 

FLAG1 = 1 ALERTS MAIN PROG. THAT THERE IS A NUMERICAL 
PROBLEM AND TO RESET H=I 

FLAG2 = 1 ALERTS SEEK THAT CHK IS RETURNING A NEW 

SEARCH VECTOR. IT IS RESET IN CHK WHEN NO 
NEW CONSTRAINTS ARE VIOLATED OR THE SEARCH 
VECTOR POINTS INTO PEAS. SPACE WITH NO 
EQUALITY CONSTRAINTS PRESENT. 

FLAG2 = 2 CAUSES FLAG4=1 AND RETURNS TO THE MAIN PROG. 

TO DO THE DFP METHOD. 

FLAG3 = 1 ALERTS PROG. THAT A BOUND EXTREMA IS LOCATED. 

FLAG3 = 2 ALERTS PROG. THAT THE LINE SEARCH CANNOT FIND A 
BETTER VALUE THAN THE LAST ONE. (MAX. STEPS) 

FLAG4 = 1 ALERTS MAIN PROG. TO INITIATE THE DFP METHOD. 

FOR SUBROUTINE CHK. 

FLAG? = 0 ALERTS CHK THAT THREE STEPS ALONG THE LINE 

SEARCH WERE MADE( MAX NO.) WITHOUT FINDING AN EXTREMA. 

FLAG? = 1 ALERTS CHK THAT AN EXTREMA WAS LOCATED ALONG 
THE SEARCH VECTOR. 

FLAGS = 1 SET WHEN GOING FROM TWO POINT FIT TO A THREE 
POINT METHOD. SET = 0 OTHERWISE. 

COMMON X0(10) ,X1(10) ,H(10, 10) ,A(20, 10) ,GR0(10) ,BC(20) , 
1N1,NC,AAT2(10,10) ,K0F,NCV0(20) ,NLEC,NEC,NGEC,NGEX,KC 
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DIl'lENSION DD(10) 

IilTUCiER PLA(;1 , P^AG2 , FojAa3 , FMG4 , FLAG6 , PLAG7 , FLAGS 
DATA ETOL1 ,ET0L2/10. , 1 .L'-5/ 

IU=0 
IC1IT=0 
?LAG2=0 
FLAG6=0 
FLAG7=0 
DO 205 1=1 , D1 
205 DD(I)=0.0 


*****DD=ii*GKAD IS TnK SEARCrf VECTOR. 


DO 150 I=1,:J1 
DO 149 J = 1,N1 

149 DD(l)=ii(l,J)*GHO(J)+DD(l) 

150 COE'i’INGE 

151 00:11’ I NUL 
FLAG3=0 

iVRITE(6,307)iM1 ,(DD(I),I=1 ,i1 ) 

307 F0RMAT(4X, 'SEARCtl VECTOR IE LliE 


o 


LARCH' /NE13.6) 


***** DETECTIOH OF TnS COHDITIOE FOR AN UNBOUE DED SOLQTIOH WILL BE 

DOEE. IF THESE CONDITION ARE MET THEN A MESSAGE WILL BE PASSED ALONG BUT 
TiiE PROG. WILL NOT BE STOPPED( COULD STILL FIND AN EXTREME VALUE) . 


ALL CONSTRAINTS IN WHICH THE 
AND NORMAL TO THE CONSTRAINT 
IDH=-1 

DO 12 r = 1 ,NC 
IF(I.GT.NLEC)IDR=1 
UN3x^lD=0.0 
DO 12 J=1 ,N1 

UNBnD=UN.3ND+DD( J)*IDR*A(I, J) 


DOT PRODUCT 0? THE SEARCH VECTOR 
PLANE ARE >/=0 ARE NOT BINDING. 


c***** at least one CONSTRAINT v/ITH A NEG. DOT PRODUCT DENOTES A BINDING 
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C CONSTRAINT HAS BEEN FOUND. 

C 

TF(UNBND.LT. 0.0)00 TO 13 
TF(I.NE.NC)GO TO 12 
WRITE(6,310) 

310 FORMATC CONDITIONS FOR AN UNBOUNDED SOLUTION HAS BEEN DETECTED') 
GO TO 13 

12 CONTINUE 

13 CONTINUE 
DMP=0.0 
DM =0.0 

IS USED TO SET THE SENSITIVITY OF THE LINE SEARCH TO THE 
MAGNITUDE OF THE HGRAD. 

DO 3 I=1,N1 
3 DMP=DMP+DD(T)»DD(I ) 

DMPrSQRT(DMP) 

PT=DMP 

DO 18 1=1 ,N1 

18 DM=DM+NUM*GRO(I)*DD(l)/DMP 

DM = THE GRAD. OF A FUNCT OF RHO. (XO+RHO*HGRAD) 

IF(DMP.GT. 10000. )RHO=. 001 

IF (DMP.GE. 200. .AND. DMP.LE. 10000. ) RHO =1 O./DMP 
IF(DMP. GT. 5. .AND.DMP.lt. 200)RH0=. 1 
IF(DMP.LE.5)RH0=1. 

WRITE (6, 300) 

300 FORMATCTHIS IS H(I,J) IN LINE SEARCH') 

WRITE(6,303)(N1 ,(H(I,J),J = 1,N1),I = 1,N1 ) 

303 F0RMAT(NE13.6) 

C 

C**»«*JF THE HESS IS BADLY CONDITIONED OR HAS EXCESSIVE ROUNDOFF ERROR 
C THE SEARCH VECTOR MAY HAVE THE WRONG SIGN. 



C HENCE DD*GRAD=+ ALWAYS UNLESS SOMETHING IS WRONG. 

this test on THE CONDITION OF SEARCHING IN THE CONSTRAINT MANIFOLD. 
C 

IF(FLAG2.EQ. 1 )GO TO 10 
IF(DM*NUM.LT.O.O)GOTO 7 
C 

STEP DOWN THE SEARCH VECTOR FOR AT LEAST THREE STEPS 
C THE LATER INTERPOLATION USES THREE POINTS 

C 

10 CONTINUE 
LLrO 

RH01=RH0 

RHOO=0.0 

RHO0O=0.0 

YO=YL 

YOO=YL 

110 CONTINUE 

DO 120 T=1,N1 
120 X1(I)=X0(I)+DD(I)*RH0 
CALL FUNCT(X1,Y1,N1) 

C 

C*»***tHE NEXT IF HAS (Y1.GE.Y0) WHEN MIN. AND (Y1.LE.Y0) FOR MAX. 

C 

IF(NUM.LT.O..AND.Y1.GE.YO)GO TO 2 
IF(NUM.GT.O. .AND.Y1.LE.Y0)G0 TO 2 
IF(LL.EQ.O)GO TO 1 
YOOzYO 
RHOOO=RHOO 
1 Y0=Y1 
RHOOrRHO 
LL=LL + 1 
C 

C»****IF MAX(MIN) HASN'T BEEN DETECTED YET THEN CALC. A PT. ANYWAY. 

C 

IF(LL.EO. 3. AND.FLAG2.E0. 1 )G0 TO 8 
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TF(LL.LE. 10)G0 TO 22 
FLAG3=2 
y=Y1 
GO TO 8 
22 CONTINUE 

RH0=RH01*2.**LL 
RHOMsRHO 
GO TO 110 

CATCHES THE CASE WERE THE FIRST PT. STRADDLES THE EXTREMA. 
2 IF(LL.E0.0)G0 TO 6 


**** THE ABOVE CATCHES CASE TWO ( y ' ( d 1 ) ,y ( cil ) ,v ( d2 ) ,d 1 , d2 ARE KNOWN) 
***** CASE ONE IS WHERE y ( d1 ) ,y (d2 ) ,y ( d3 ) , d1 , d2 , d3 ARE KNOWN. 


19 CONTINUE 

XX=(RHOO-RHOOO)*(Y1-YOO)~(RHO-RHOOO)*(YO-YOO) 
RHOM=RHOOO+.5*( ( ( RHOO-RHOOO ) **2* (Y 1 -YOO )- (RHO-RHOOO ) **2 * 
X(YO~YOO))/XX) 

IF(ABS(XX).LT. .0000000001 )RHOM=. 000001 

4 CONTINUE 

DO 5 1=1, N1 

5 X1(I)=X0(I)+DD(I)*ABS(RH0M) 

WRITE(6,308) 

308 F0RMAT(4X, 'XI(I)' ,9X,'X0(T)' ,9X,'DD(T)' ) 

WRITE (6, 306) (XI (I ) ,X0 (I ) , DD (I) ) , T=1 ,N1 
306 F0RMAT(3E13.6) 

CALL FUNCT(X1,Y,N1 ) 

FLAG7=1 alerts CHK THAT AN EXTREMA ON THE LINE SEARCH 
WAS FOUND 


VO 



PLAa7=1 

IP(NUM.L'J?.O.AND.Y.LE.YO)GO TO 8 
IP(NUM.GT.O.AWD.Y.(}S.YO)GO TO 8 
C 

C***** GIVE IT TEN TRIES AND THEN CNEGK TO SEE IP ITS CLOSE ENOUGH. 
C 

ICNT=ICNT+1 

IP(ICNT.GT.10)G0 TO 15 

IP(PLAG8.EQ.1 )GO TO 17 

IP(LL.NE.O)GO TO 14 

RriOO=RiIOM 

YO=Y 

PLAG8=1 

GO TO 19 

14 IP(ABS(RHOM) .GT.RHOO)GO TO 20 
Rri000=RH0M 

YOO=Y 
GO TO 19 
20 CONTINUE 
RriO=RriOM 
Y1 =Y 

GO TO 19 
17 RxlO=RHOO 
Y1=Y0 
RH00=RH0M 
Y0=Y 

GO TO 19 

15 CONTINUE 

IP(ABS(ABS(Y)-ABS(YO) ) .LT.1 .E-5*ABS(Y0) )G0 TO 8 

16 CONTINUE 
PLAG3=2 
WRITE (6, 31 2) 

312 FORMAT ( 'SEARCH FAILED TO FIND A BETTER VALUE THAN XO') 

CALL GRAD(X0,GR0,N1 ,NUM) 

GO TO 9 
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***** THIS JS THE FINAL ADJUSTMENT OF THE SEARCH STEP SIZE 

6 TF(ABS(Y1).LT.(100.*ABS(YO)))30 TO 21 
RHOr. 1 *RH0 

GO TO 10 
21 CONTINUE 

RHOM = .5*RHO**2*DM/(RHO*DiM+YO-Y1 ) 

IF( ABS(RH0*DM+Y0-Y1 ) .LT. .0000000001 )RH0M=. 000001 

***** jF the grad or PROJECTION ARE LARGE AND THE SEARCH STEP SIZE IS 
SMALL THEN MORE ACCURACY IS CALLED FOR THAN THE QUADRATIC FIT 
CAM GIVE SO SUBROUTINE FIB IS CALLED. 

IF(PT.GT.ET0L1.AND.RH0M.LT.ET0L2)CALL FIB( RHO, RHOM , NUM , XO , N1 , DD ) 
GO TO 4 

7 FLAG1=1 

8 CONTINUE 
IF(NC.EQ.O)GO TO 9 

***** this ROUTINE CHECKS TO SEE IF ANY CONSTRAINTS ARE VIOLATED 
AND IF SO IT RETURNS WITH A NEW SEARCH VECTOR 

CALL CHK (PT , DD , FLAG2, FLAGS, FLAG6, NUM , FLAG?, 10 ) 

***** fLAG3=1 A CONSTRAINED MAX(MTN) HAS BEEN FOUND 

IF(FLAG3.EQ. 1 )G0 TO 9 

***** fLAG2=1 CONSTRAINT(S) have been VIOLATED AND CHK HAS 
SUPPLIED A NEW SEARCH VECTOR. 

FLAG7=0 

IF(FLAG2.NE. 1 )GO TO 11 

CALL FUNCT(X0,YL,N1 ) 


Os 



ICNTrO 
GO TO 151 

11 IF(FLAG2.EQ.2)FLAGi» = 1 
9 CONTINUE 
RETURN 
END 
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SUBROUTINE FTB( RHO , RHOM , NUM , XO , N1 , DD ) 

THIS IS AN 8 STEP EIBONACCI SEARCH 

DIMENSION XO(N1) ,X1 (10) ,XT1(10) ,DD(N1 ) 

DATA ET0L/1.E-5/ 

REAL L2 
ISTEPrI 
RANU=RHO 
RANL=0. 

L2=(RH0*13.+(ET0L*(-1 )»*8))/21. 

X2=RH0-L2 

1 CONTINUE 

IF(ISTEP.EQ.8)G0 TO 13 

ISTEP=ISTEP+1 

DO 2 1=1, N1 

XT1 (I )=XO(I )+DD(I)*L2 

2 X1(I)=X0(I)+DD(I)*X2 
CALL FUNCT(XT1 ,Y1,N1 ) 

CALL FUNCT(X1 ,Y2,N1 ) 

IF(Y1.LT.Y2.AND.X2.GT.((RANU-RANL)/2.+RANL))G0 TO 10 
IF(Y1.LT.Y2.AND.X2.LT.( (RANU-RANL)/2.+RANL) )G0 TO 7 
IF(Y1.GT.Y2.AND.-X2.GT.((RANU-RANL)/2.+RANL))G0 TO 6 
IF(NUM.GT.O)GO TO 12 

3 RANU=L2 

4 CONTINUE 

IF(X2.GT.((RANU-RANL)/2.+RANL))G0 TO 5 

L2=(RANU-X2)+RANL 

GO TO 1 

5 L2=RANU-(X2-RANL) 

GO TO 1 

6 CONTINUE 

IF(NUM.GT.O)GO TO 11 
14 RANL=L2 
GO TO 4 


0 \ 


0 \ 


7 CONTINUE 
JF(NUM.GT.O)GO TO 3 

12 RANL=X2 

8 CONTINUE 

IE(L2.GT.((RANU-RANL)/2.+RANL))G0 TO 9 

X2=(RANU-L2)+RANL 

GO TO 1 

9 X2=RANU~(L2-RANL) 

GO TO 1 

10 CONTINUE 
IF(NUM.GT.O)GO TO U 

11 RANU=X2 
GO TO 8 

1 3 CONTINUE 

RH0Md(RANU-RANL)/2.+RANL 

RETURN 

END 
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SUBROai'Ii'lE CrilC( P'i’ , P , PLAG2 , PLAG3 , FLAG6 , EW , FijAG7 , 10 ) 
G*****im - THE X^IUHBER 0? TRE VIOLATED GONSTRAIRT 

DCS - THE NUMBER OP TxlE VIOMTED CONSTRAINTS (NCV(l)) 
TxiAT RAVE TRE LOWEST ALPHA. THIS IS USED TO 
INDEX NCV. 

***** j3X: 

NCV(NCS(2)) - IS THE SECOND OP THE NEWLY VIOLATED 
CONSTRAINTS THAT RAVE ALPRA'S EQUAL TO THE LOWEST 
ALPHA (MINIMUM SET VIOLATED). 

PLAG2 = 0 ALERTS CHK THAT TRE PREVIOUS CYCLE DIDN'T 
VIOLATE A CONSTRAINT. 


PLAG2 


ALERTS CHX TiIAT CONSTRAINTS RAVE ALREADY BEEN VIOLATED. 


PLAG3 = 1 ALERTS PROG. THAT A BOUND EXTREMA HAS BEEN POUND. 

PLAG6 > 0 ALERTS CHK THAT DROP tiAS GENERATED A NEW POINT. 

PLAG6 = 2 ALERTS CHK TO LOOK FOR THE PURTREREST CONSTRAINT 
PROM POINT RETURNED BY DROP AS IT IS IN 
INFEASIBLE SPACE (SEARCH VECTOR POINTS IN). 

FLAG? = 1 ALERTS CHK THAT AN EXTREMA WAS POUND ON LINE SEARCH. 


COMMON X0(10),X1 (10), H( 10,10) ,A( 20,1 0),GR0(10),BC(20) , 
INI, NC , AAT2 ( 1 0 , 1 0 ) , KOP , xJC VO ( 20 ) , NLEC , NEC , NGEC , NGEX , KC 
D IMENSION B ( 20 ) , NG V ( 20 ) , SAX ( 20 ) , SAxiD ( 20 ) , ALPHA ( 20 ) , 
1NGS(20),XN(10),P(10),GRN(10),HGR0(20),NT0C(20), 

2MNB ( 20 ) , NC VI ( 20 ) , GRT ( 1 0 ) , XT ( 1 0 ) 

I NTEGER FLAG1 , PLAG2 , PLAG3 , PLAG4 , PLAGo , FLAG? 

REAL IDENT 
DATA ET0L/1.E-5/ 

WRITE(6,313) 

313 FORMAT ('THE XI POINT') 
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WRITE(6, 31^)N1 , (X1(T ) ,T=1 ,N1 ) 

F0RMAT(NE1 3.6) 

KF=0 
FLAGi|=0 
BT=0. 0 
SAX(0)=0.0 
SAHD(0)=0.0 
DO U I=1,NC 
B(I ) = 0.0 
HGROd )=0.0 
1i| CONTINUE 

***** THIS CALCULATES THE ACTUAL VALUES OF THE CONSTRAINT EQ.'S 

DO 1 I=1,NC 
DO 1 J=1,N1 

B(I)=A(I,J)*X1(J)+B(I) 

IF(I.GT.N1)G0 TO 1 
IF(FLAG2.EQ. 1)G0 TO 1 

*#*«*0NLY NEED HGRO WHEN GOING FROM FEASIBLE SPACE TO CONSTRAINT. 


HGROd ) = H(I, J)*GR0(J)+HGR0(I ) 
1 CONTINUE 


SORT THE CONSTRAINTS FOR THOSE THAT ARE VIOLATED. 

NEWLY VIOLATED CONSTRAINTS 
OLD CONSTRAINTS STILL VIOLATED 
CONSTRAINT VIOLATIONS 
THE OLD VIOLATED CONSTRAINT 
" " " '• STILL VIOLATED 

" NEWLY VIOLATED CONSTRAINT 

****»*****K***¥t*it**itjlc1iit*****iiie*»1HiiK**1t******1c****iHi*it**iHtiHt** 


KA 

= 

NO. 

OF 

KO 

= 

n 

II 

KOF 

= 

II 

II 

NCVOd 

) = 

II 

It 

NCV1 (I 

) = 

II 

II 

NCVd ) 

= 

It 

If 


Q*******ikic*iiii**iii*ii**iiitit**it****¥iiHc*iHt*itiie***itiiii********ie*itik***1iii 
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KA=0 

K0=0 

DO 2 T=1,NC 
BF=BC(I)-B(I) 

C 

this is here to alleviate numerical difficulties 

c 

IF(ABS(BF) .LT.ETOL.AND.FLAG2.EO. 1 )GO TO 51 
IF(FLAG6.NE.2)GO TO 62 

QiHH»*«*»**»***#**»*Ji**fe******fc******'********************'********** 

c* ***************************************************************** 

THE FOLLOWING SECTION IS HERE TO CHEK FOR A NEW CONSTRAINT 
BEING ENCOUNTERED WHEN THE SEARCH VECTOR POINTS INTO 
FEASIBLE SPACE AFTER THE PROGRAM HAS STEPPED INTO 
INFEASIBLE SPACE. 

***************************************************************** 

fcjk***»»««*»*»***************if*******»‘******'**S'***********^******** 

IF(BF.LT.O. .AND.I.LE.NLEOGO TO 2 
IF(BF.GT.O. .AND.I.GT.NLEC+NEOGO TO 2 
DO 6M IIr1,K0F0 
WRITE(6,328)I,NT0C(II) 

328 F0RMAT('I=' ,13,' NT0C=',I3) 

64 IF(I.EQ.NT0C(II))G0 TO 46 
GO TO 2 
62 CONTINUE 

IF(BF.GE. 0.0. AND.I.LE.NLEOGO TO 2 
IF(BF.LE.O.O.AND.I.GT.(NLEC+NEC))GO TO 2 
51 CONTINUE 

IF(KOF.EQ.O)GO TO 46 
DO 45 K=1,KOF 

45 IF(NCVO(K).EO.I)GO TO 50 

46 CONTINUE 
KA=KA+1 
NCV(KA)=I 
GO TO 2 


OS 
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0 \ 

00 


50 K0=K0+1 
NCV1(K0)=T, 

2 CONTINUE 

***** THIS TESTS TO SEE IF ANY NSW CONSTRAINTS WERE VIOLATED 
AFTER THE BRANCH IT TESTS TO SEE IF IT WAS ALREADY 
IN A CONSTRAINT PLANE(SEE COMMENT BELOW). 

IF(KA.EQ.O)GO TO 13 

WRITE ( 6,300 )KA, (NCV(I) ,1=1 ,KA ) 

300 F0RMAT(I3,' CONSTRAINTS HAVE BEEN VIOLATED. THEY ARE', 2014) 

***** FLAG4 IS SET WHEN THE NEXT XO BEING LOOKED FOR IS DOWN A 
SEARCH VECTOR ALREADY IN A PLANE THAT HAS VIOLATED AT 
LEAST ONE OTHER PLANE. 

***** prom HERE TO STATEMENT 4 A SIMULTANEOUS SOLUTION OF 

THE SEARCH VECTOR AND THE CONSTRAINT PLANES IS CALCULATED 
THIS ACCOMPLISHES TWO THINGS 

(1) AN EXACT PT. (XI) IN THE PLANE IS LOCATED 

(2) AN EXACT DISTANCE ( ALPHA ) FROM XO TO XI IS CALCULATED 

*********«****************1c************************************* 

FLAG4=FLAG2 
FLAG2=1 
FLAG3=0 
DO 15 1=1 ,KA 
SAX(I )=0. 0 
SAHDd ) = 0.0 
15 CONTINUE 
DO 4 1=1 ,KA 
DO 3 K=1,N1 

SAXd ) = SAX(I)+A(NCV(I ) ,K)*XO(K) 

IF(FLAG4.E0.0 )G0 TO 24 

SAHDd ) = SAHD(I )+N UM* A ( NC V (I ) , K ) *P (K ) 
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GO TO 5 
24 CONTINUE 

SAHD ( I ) =SAHD ( I ) +NJM*A ( i'fC V ( I ) , K ) *HGR0 ( K ) 

3 CONTINUE 

4 ALPHA(I)=(BC(NCV(I))-3AX(I) )/SAHD(l) 

C 

C*****TiiIS ALPHA IS IN THE EQ. XI =XO+ALPHA*HGRO .IT WILL BE USED 
TO ASCERTAIN THE SMALLEST BORDER OP THli FEASIBLE SPACE 
THIS WILL BE DONE BY SORTING FOR THE SMALLEST ALPHA 

DO 16 1=1 ,KA 
DO 16 K=1 ,N1 

****»SORT FOR THE LOWEST ALPHA OF THE NEWLY VIOLATED CONSTRAINTS 
***** XP PLAG6 = 2 FIND THE LARGEST ALPHA ( SEARCH PROM INPEAS. SPACE) 

TRIAL=1 .E-7 

IP(PLAG6.GT.O)TRIAL=1 .E-14 
NM = 1 

NC3(1 )=1 

ALPO=ABS( ALPHA (1 ) ) 

DO 5 1=1 ,KA 

IF(PLAG6.GT.O)GO TO 38 

IF(ABS(ABS(ALPHA(I) )-ALPO) .LT.TRIADGO TO 52 

38 CONTINUE 

IP(PLAG6.NE.2)GO TO 43 
IF(ABS(ALPHA(l)).LT.ALPO)GO TO 5 
GO TO 52 
43 CONTINUE 

IF(ABS(ALPHA(I) ) .GT.ALPO)GO TO 5 
52 IF(I.EQ.1 )GO TO 5 
ALPO=ABS(ALPHA(l) ) 

NM=I 

5 CONTINUE 
K=1 



IF( PLAG6 . EQ . 2 ) FLAa6=0 
NPLAG1 0=0 
C 

G*****0NLY ONE CONSTRAINT NAS THE LOWEST ALPHA . 

C 

IP(NM.EQ. 1 )G0 TO 7 
C 

C*****HOW MANY OE THESE CONSTRAINTS HAVE ALPHA'S = TO THE LOWEST 
C ALPHA AND WHA‘ IS THEIR NUMBER ( SET NCS) . 

C 

DO 6 1=1, NM 

IE(ABS(ABS(ALPHA(I) )-ALPO) .GT.TRIAL)GO TO 6 

NCS(K)=I 

K=K+1 

6 CONTINUE 
IP(NM.EQ.2)NCS(2)=2 
K=K— 1 

7 CONTINUE 
KP=K 

Q *********** ************** ****************************** ****** 

Q************** ********************************** ************** 

c 

C THIS SECTION EXAMINES THE NEWLY ACQUIRED CONSTRAINTS. 

C IP IT DETERMINES THAT IT HAS POUND A NEW ONE THEN 

C PLAG6 IS RESET TO 0. PLAG6 IS NOT ZERO IP IT HAS BEEN 

C THROUGH DROP. NTOC CONTAINS ALL OP THE PREVIOUSLY 

C VIOLATED CONSTRAINTS. 

C 

Q** ************************************************ ************* 
Q ******************************************************* ******** 

IP(PLAG6.EQ.O)GO TO 63 
DO 68 1=1 ,KP 
DO 66 11=1, KOPO 

66 IP(NCV(NCS(I)) .EQ.NTOC(II) )G0 TO 68 
WRITE(6,335)NGV(NCS(I) ) 
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335 F0RMAT(I3,' DOES NOT MATCH THOSE THAT WERE DROPPED') 
WRITE(6,381 )K0F0, (NTOC(II) ,TI. = 1 ,K0F0) 

FLAG6=0 
GO TO 63 

68 CONTINUE 
IF(KO.EQ.O)GO TO 73 
DO 69 T = 1,K0 

DO 67 TI=1,K0F0 

67 IF(NCV1(I).EO.NTOC(TI))GO TO 69 
WRTTE(6, 335)NCV1 (T ) 

WRITE (6, 381 )KOFO, (NTOC(II) ,11=1 ,KOFO) 

FLAG6=0 
GO TO 63 

69 CONTINUE 


IF ALL OF THE SAME CONSTRAINTS ARE RE-ACQUIRED 
THEN ASSUME A CONSTRAINED EXTREMA AT XO 

kxkkkkkkkirkkkxirkkkkkkkkkkkkkxxkxkkxkxxxkkkxkkkxxkxkkkxkkkirkkkk 

73 CONTINUE 

IF(KF+KO.NS.KOFO)GO TO 63 
CALL FUNCT(XT,I0,N1,NUM) 

CALL GRAD(XT,GRO,N1 ) 

FLAG3=1 
RETURN 
63 CONTINUE 

iiniiik ONLY THOSE CONSTRAINTS WITH ALPHA'S EQUAL TO THE 
LOWEST ALPHA WILL BE ADDED TO THE MANIFOLD. 

WRITE (6, 304 )ALPO, K, (NCV(NCS(I ) ) ,1=1 , K) 

304 FORMATCTHE SMALLEST ALPHA=' ,E1 3. 6, 

1' FOR ',13,' CONSTRAINTS. ' /'THEY ARE', 2013) 
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C***** I'llESE XO'S ARE M THE CONSTRAMi PLANE. 

C 

DO 8 1=1 ,N1 

IP(PLA04.EQ.0)00 TO 25 
XO(I)=XO(l)+ABS(ALPHA(NM) )*P(I) 

GO TO 8 
25 CONTItiUE 

XO ( I ) =XO ( I ) +AB3 ( ALPRA ( NM ) ) *RGRO ( I ) 

8 CONTINUE 

WRITE(6,410)N1 ,(X0(I),I=1 ,N1 ) 

410 PORMAT('TxRE NEW XO CALC. WITH ALPA IS ' /^'fEI 3 • 6 ) 
C ALL GRAD ( XO , GRO , N 1 , N UM ) 

CALL PUNCT(X0,Y0,N1 ) 

WRITE(6,303)YO,N1 , (GxRO(l) ,1=1 ,N1 ) 

305 PORMAT( 'P(XO)=' ,E13.6/'GRAD P(X0)= ' ,NE1 3. 6// ) 
18 CONTINUE 
KA=0 
KS=0 
KSS=0 

IE(KOF.NE.O)GO TO 47 

K0P=1 

MNB(0)=0 

IP(KP.GE.N1 )G0 TO 9 
GO TO 49 
47 CONTINUE 


***** SCAN THE DIFFERENCE BETWEEN THE OLD VIOLATED CONSTRAINTS 


AND THOSE THAT ARE STILL BEING VIOLATED 


DO 27 1=1, KOF 
DO 26 J=1 ,K0 

26 IP(NCV0(I) .EQ.NCV1 (J))G0 TO 27 
KS=KS+1 
C 

C***** KSS COUNTS ONLY THE NO. OP OLD CONSTRAINTS TO BE 
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C DROPPED THAT ARE DETECTED BY THE COHSTRAIiT TEST 

C 

KS3=KSS+1 
i'U>iB(KS) = I 
27 COWTIHaE 
C 

0 ***** EXCESS COHSTRAIHTS WILL BE DETECTED AMD DROPPED 


** tt*********i^***ilt********1Ht************1h********************* 
****************************** )^***** ************************* 

HERE WE ARE COMCERMED THAT SPACE IS COMPx,ETELY 
SPAMMED AMD THE PROJECTIOM VECTOR IS TOTALLY DEFINED 
LEAVING MO ORTHOGONAL COMPONENT. DROP WILL DETERMINE 

THOSE CONSTRAINTS TO BE DROPPED. 

******** *■*■*•******■*•*• »■*■***•****♦** ***■)(■**»*****•»•*******■*•****•»•***■*■* 

■JHt*****-********************** ******************************** 

IF( (KF+KO-KS) .LT.M1 )GO TO 55 
9 CONTINUE 

CALL DROP ( NO , NC V 1 , KF , MC V , NC S , FLAG3 , FLAG6 , P , XT , NT OC , KOFO , N UM ) 
WRITE(6,381 )KOFO, (MTOC(I) ,1=1 ,X0F0) 

381 FORMAT ( 'IP' ,NI3, ' ARE REAQUIRED AN EXTREMA IS ASSUMED') 
IP(PLAG3.EQ.1 )G0 TO 65 
IF(PLAG6.EQ.3)GO TO 21 

*************************************************************** 

* 


IF PLAG6=3 THEN THE SEARCH VECTOR POINTS INTO FEASIBLE 
SPACE. IF THERE ARE EQUALITY CONSTRAINTS PRESENT THEM 
THEY ARE KEPT AND PROJECTIONS ARE MADE ON THEM. IF 


NONE ARE PRESENT TilEN RETURN TO THE DFP METrlOD IS MADE. 


* 

* 

* 

* 

* 


*************************************************************** 


RETURN 
55 CONTINUE 

IF(NEC.EQ.O.OR.KS.EQ.O)GO TO 71 
L=0 


CO 


oooo ooooo 


'J 


KIND1=1 

WE MUST SCAN THE LIST OF CONSTRAINTS TO BE DROPPED TO DETERMINE 
IF ANY OF THEM ARE EQUALITY CONSTRAINTS. IF ONE IS THEN 
IT WILL DROPPED FROM THE LIST. 

54 CONTINUE 

DO 70 I=KIND1,KS 
DO 70 J=1,NEC 

IF(NCVO(MNB(T+L ) ) .NE. (NLSC+J ) )GO TO 70 
IF(I.EQ.KS)GO TO 59 
L=L + 1 

DO 61 KIND=T,KS~1 
61 MNB(KIND)=MNB(KIND+1 ) 

GO TO 58 

70 CONTINUE 
GO TO 71 

58 KS=KS-1 
KIMDIrl 
GO TO 54 

59 KS=KS~1 

71 CONTINUE 
IF(KS.EQ. 0)GO TO 49 
IF(KO.GT.O)GO TO 12 

WRITE (6, 321 )KS, (NCV(NCS(T )) ,1=1 ,KS) 

GO TO 49 

72 CONTINUE 

WRITE (6, 321 )KS,(NCVO(MNB(I)),I=1 ,KS) 

321 FORMAT(I3,' CONSTRAINTS WILL BE DROPPED. THEY ARE ',2014) 

49 CONTINUE 

IF(KO.EQ.O)GO TO 29 

**^********it*it1i***k**it**itiik***ik****^*ic*¥t*ii**ile*ikiiik****it*ik*****^*ii^^* 

AT THIS POINT WE UPDATE THE NCVO ARRAY BY 
DROPPING OFF THE CONSTRAINT NOT STILL 
VIOLATED. SINCE WE ASSUME ALL PLANES 
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TO BE LINEARLY INDEPENDENT THEN THE NUMBER OF 
CONSTRAINTS VIOLATED CAN NEVER BE GREATER THAN 
THE DIMENSIONS OF SPACE. 

iH*it»*»**»iiiiim»**ii»*»**:**ii»ii***it***i>tk*ili»iHi*ii*it*7i******i^******** 

IDDrO 

DO 28 1=1 ,KO 
IF(KS.E0.0)G0 TO 57 
IF(KSS.EO.O)KSS=1 
DO 56 J=KSS,KS 

56 IF(NCV1(I).E0.NCVO(MNB(J)))GO TO 28 

57 CONTINUE 
IDD=IDD+1 
NCVO(IDD)=NCV1(I ) 

28 CONTINUE 
KO=KO-KS 

29 CONTINUE 
IF(KS.EQ.0)G0 TO M8 
NSN = -1 


***** this ROUTINE CALCULATES AND UPDATES THE PROJECTION MATRIX(AAT2) 
THE ARGUMENT MNB HAS NO EFFECT IN SPAN FOR SUBTRACTION 

CALL SPAN(NCVO,MNB,IDD,NSN,IO) 
il8 CONTINUE 

IF(KF.EQ.O)GO TO 39 
DO 44 1=1, KF 

44 NCVO(KO+I ) = NCV(NCS(I )) 

WRITE (6, 322 )KF , (NCV (NCS (I ) ) , I = 1 , KF ) 

322 FORMAT(I3,' CONSTRAINTS WILL BE ADDED. THEY ARE', 2014) 

NSN=1 

CALL SPAN(NCV,NCS,KF,NSN,IO) 

39 CONTINUE 
KC=KC+1 
DO 31 I=1,N1 
31 P(I)=0.0 


Ui 


oooo oooo 


WRITE (6, 333) 

333 FORMATCTHIS IS THE PRESENT PROJECTION MATRIX') 
WRITE(6, 314) (N1 ,(AAT2(I, J) , J=1 ,N1) ,1=1 ,N1 ) 

DO 40 J=1 ,M1 
DO 40 I = 1 , N 1 

***** P IS THE PROJECTION OF THE GRAD ON THE CONSTRAINT 
MANIFOLD. IT IS THE NEW SEARCH VECTOR 

IF(I.NE.J)GO TO 16 

P(I )=(1.-AAT2(I,J))*GR0(J)+P(I ) 

GO TO 40 

16 P(I)=-AAT2(I, J)*GR0(J)+P(I) 

4 0 CONTINUE 
PT=0.0 
DO 42 I=1,N1 
42 PT=PT+P(I )**2 
KOF=KF+KO 
PT=SQRT(PT) 

***** IF PT=0 AND GRAD DOES NOT POINT INTO FEASIBLE SPACE 
THEN A CONSTRAINED MAX(MIN) HAS BEEN MET. 

IF(PT.GT. .001 )GO TO 19 
FLAG3=1 
GO TO 65 
19 CONTINUE 

Q**********it****¥titit*ikiit*ic**iiitititic***i(*1tii*k***^¥i*****it***ic***** 


C * 
C IF WE HAVE DONE A SMALL STEP INTO FEASIBLE SPACE * 
C THEN WE WILL CHECK AT THIS POINT TO DETERMINE IF * 
C THE NEW CONSTRAINT MANTFOLLD CONTAINS THE * 
C EQUALITY CONSTRAINTS. THIS NEED ONLY BE DONE AT * 
C A POINT WHERE THE PROJECTION VECTOR POINTS * 
C DIRECTLY AWAY FROM THE ORIGINAL POINT WHERE ALL * 
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OF THE CONSTRAINT PLANES INTERSECTED. THIS IS * 

DETERMINED BY: * 

+1 =(XO-XT)/ABS(XO-XT)*P/PT * 

IF(FLAG6.EQ.O)GO TO 11 
XTDzO. 

XTDTrO. 

DO 12 Ir1,N1 

12 XTD=XTD+(XO(I )-XT(I) )**2 
XTDrSQRT(XTD) 

DO 30 1=1, N1 

30 XTDT = (X0(I )-XT(I))*P(I )/(XTD*PT)+XTDT 

***** CHECK TO SEE IF ALL CONSTRAINTS ARE GENERATED 

IF(XTDT.LT.0.0)G0 TO 11 
IF(ABS(ABS(XTDT)-1.).GT. 1.E-12)G0 TO 11 
WRITE(6, 387)ABS(ABS(XTDT)-1 . ) 

387 FORMATCMUST BE > THAN 1.E-12 BUT IS',EU.6) 

***** SCAN MANIFOLD TO SEE IF ALL EQUALITY CONSTRAINTS ARE PRESENT 

M=0 

IF(KF.EQ.0)G0 TO 33 
DO 32 I=1,KF 

IF(NCV(NCS(I) ).LE.NLEC.OR.NCV(NCS(D) ,GT. (NLEC+NEC))GO TO 32 
M=M + 1 

32 CONTINUE 

33 IF(KO.EQ.O)GO TO 37 
DO 34 I=1,KO 

IF(NCV1(I).LE.NLEC.OR.NCV1(I).GT.(NLEC+NEC))GO TO 34 
M=M+1 

34 CONTINUE 
37 CONTINUE 
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00 


c***** all EQUALI'j?!f CONSTRAIUl’S ARE PRESENT 'rfilEN H=NEC 
C 

IP(M.£Q.KEC)GO TO 36 
DO 35 1=1, ill 

35 X1(I)=XT(I) 

PLAG3=1 
RETURN 

36 PLAG6=0 
KUE0=0 

11 CONTINUE 
RETURN 

13 IF(ELAG2.EQ.1 )G0 TO 17 
WRITE (6, 302) 

302 ?ORMAT('NO CONSTRAINTS NAVE BEEN VIOLATED') 

10 PLAG2=0 
RETURN 
17 CONTINUE 

***•*•****• K-******** ************************************* )(■****** 

KU=0 TRE CASE WHERE THE GRAD POINTS INTO PEAS. SPACE * 

* 

KU=KOP NO CONSTRAINTS WILL BE DROPPED PROM AAT2 * 

★ 

KO<KOP SOME CONSTRAINTS ARE TO BE DROPPED * 

* 

PLAG7=1 AN EXTREMA WAS POUND ON THE LINE SEARCH * 

**************************************************************** 

IP(KO.EQ.O)GO TO 21 
CALL GRAD(X1 ,GRT,N1 ,NUM) 

IP(PLAG7.EQ.1 )G0 TO 74 
IP(KO.LT.KOP)GO TO 18 
DO 79 1=1 ,N1 
79 X0(I)=X1 (I) 

RETURN 

74 DO 78 1=1, N1 
78 GR0(I)=GRT(I) 


If (KO.EQ.O)GO TO 29 
DO 27 T=1,KO 
XI (T) = 0. 

AMT(T ) = AMV 

TF(NCV1(T ).LE.NLEC+NEC)AMT(T ) = -AMV 
DO 24 Jr1,N1 

24 X1(T)=X1(I)+A(NCV1(T),J)**2 
27 XI (I)=SQRT(X1 (I)) 

DO 25 1=1 ,KO 
DO 25 J=1,N1 

X0(J) = X0(J)+A(MCV1 (I) ,J)/(X1 (T )^-AMT(T )i-5.E5) 

25 CONTINUE 
29 CONTINUE 

If (Kf .EO.O)GO TO 48 
DO 33 I=1,KF 
X1(.T. )=0. 

AMT(I )=AMV 

If (NCV(NCS(I )) .LE.NLEC+NEOAMTd ) = ~AMV 
DO 32 J=1,N1 

32 XI (I )=X1 (I )+A(NCV(NCS(I ) ) , J'd*2 

33 X1(I)=SQRT(X1(I)) 

DO 34 I=1,Kf 

DO 34 J=1,N1 

XO ( J ) = X0 ( J )+A ( NC V (NCS ( I ) ) , J ) /( X 1 ( I ) '-AMTC I ) *5 . E5 ) 

34 CONTINUE 
48 CONTINUE 

DO 35 1=1, N1 
NCVO(I)=0.0 

35 NCV1 (I )=0.0 

NRITE(6,301 )M1,(XO(I),I = 1 ,N1 ) 

301 FORMATCDE WILL START AT THE FOLLOWING POINT & FOLLOW THE GRAD' 
*/NE1 3.6) 

RETURN 

END 


'4 

VO 
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00 

o 


IP(KO.LT.KOP)aO TO 18 
DO 23 1=1, N1 
23 X0(I)=X1 (I) 

GO TO 39 

21 IP(iJiSC.i\[E.O)GO TO 75 
PLAG2=2 
X0P=0 

65 CONTINUE 
DO 53 1=1, I'll 
XI (I)=XO(I) 

DO 53 J=1,N1 
53 AAT2(I, J)=0.0 
RETURN 

***-******************************-)(-*^HHt*** ****** *********** 


IP THE GRAD POINTS INTO FEASIBLE SPACE IT NILL * 
STILL HAVE TO SATISFY THE EQUALITY CONSTRAINTS * 
HENCE THE GRAD WILL CONTINUE TO BE PROJECTED * 
ONTO THE EQUALITY CONSTRAINTS. ONLY THE EQUALITY * 
CONSTRAINTS r'lANIPOLD WILL BE KEPT. * 


C***-**»*-***************)t**** ************************** ********* 

75 DO 76 1=1 ,x'i1 
DO 76 J=1 ,N1 

76 AAT2(I,J)=0.0 
10=0 

K0=0 

KP=NEC 

DO 77 1=1 ,KP 
iNCV(I)=xNLEC+I 
NCS(I)=I 
NCV0(I)=NLEC+I 

77 NCV1 (I)=NLEC+I 
GO TO 48 

END 
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SU3R0UT JWE SPAN (NO , NOC , KN , NSN , TO ) 

THIS ROUTINE CALCULATES AND UPDATES THE PROJECTION HATHIX(AAT2) 
IT USES THE GRAM-SCHMIDT METHOD AND IS WRITTEN IN THE 
USUAL MATRIX NOTATION. 

COMMON X0(10) ,X1 (10) ,H(10, 10) ,A(2 0, 10) ,GR0(10) ,BC(20) , 

1 N 1 , NC , A AT2 ( 1 0 , 1 0 ) , KOE , NC VO ( 20 ) , NLEC , NEC , NGEC , NGE X , KC 
DIMENSIO N AL ( 1 0 , 1 0 ) , NOC ( 1 0 ) , NO ( 1 0 ) , I D PAST (10) 

REAL I DENT 
C 

Ck»ii*x Xf NSW<0 SET AAT2=0.0 THEN RECALCULATE AAT2 WITHOUT 
C THE DROPPED CONSTRAINTS. 

lEdO.EO. 0)11=0 
10 = 1 
JK = 1 

IF(NSN.LT.0)G0 TO 10 
12 CONTINUE 
DO 1 J=JK,KN 
DO 1 X=1,W1 
1 AL(J,K)=0-0 
9 CONTINUE 
DO 7 I=JK,KN 
IF(N3N.GT.O)GO TO 14 
NOC(I)=I 
14 CONTINUE 
11=11+1 

IDPASTdl )=NO(NOC(I ) ) 

W R I T E ( 6 , 3 0 4 ) N 1 , ( A ( N 0 ( N 0 C ( I ) ) , J ) , J = d M 1 ) 

304 FORMATCTHIS IS THE PRESENT CONSTRAINT EO. IN SPAN ' /NE 1 3. 7 ) 

DO 3 J=1,N1 
DO 3 K=1 ,M1 
IF(J.NE.K)GO TO 23 


00 



AL(J, J) = (1 .0-AAT2(J,K) )*A(NO(NOC(T)) ,K)+AL(T, J) 

■ GO TO 3 

23 AL(T, J)=-AAT2(J,K)*A(MO(MOC(T )) ,K)+AL(T, J) 

3 CONTINUE 
ALABS=0.0 
DO 4 K=1,N1 

4 ALABS = ALABS+AL(I,K)*’'2 
ALABS=SQRT( ALABS) 

IF(ALABS.LT. .0000000001 )GO TO 7 
DO 5 J=1,N1 

5 AL(I, J)aAL(I, J)/ALABS 
DO 6 K=1,N1 

DO 6 L = 1 ,I'11 
TF(L.LT.K)GO TO 19 

AAT2(K,L)=AAT2(K,L)+AL(T,K)*AL(I,L) 

19 AAT2(L,K)=AAT2(K,L) 

6 CONTINUE 

7 CONTINUE 
RETURN 

SUBTRACT DIMENSIONS TO LOWEST ONE TO BE REMOVED 

10 CONTINUE 

DO 11 J = 1 , N 1 
DO 11 K = 1 , N 1 

11 AAT2(J,K)=0.0 
DO 2 J=1,KN 

2 IF(NO(MOC(J)).LT.(II-KN))GO TO 16 
DO 17 I=TI,II~KN,~1 
DO 17 K = 1 , N 1 
DO 17 L=1,M1 
TF(L.LT.K)GO TO 20 
AAT2(K,L)=AAT2(K,L)-AL(T,K)«-AL(I,L) 

20 AAT2(L,K)=AAT2(K,L) 

17 CONTINUE 

RETURN 
16 CONTINUE 



o o o o 


JK = I.T 

DO 8 J=1 , JJ 
DO 8 K=1,KN 

IF(NO(NOC(K)).NE.IDPAST(J))GO TO 8 

IF( J.GE. JK)GO TO 8 

JK=J 

8 CONTINUE 

JiC IS EQUAL TO THE LOWEST AL TO BE DROPPED. THE ONES BELOW THIS 
NEED NOT BE RECALCULATED. 

IF(JK.LT.2)G0 TO 15 

DO 13 I = 1,JK-1 

DO 13 K=1 ,N1 

DO 13 L=1,N1 

IF(L.LT.K)GO TO 21 

A AT2 ( K , L ) =A AT2 (K , L ) +AL ( I , K ) * AL ( I , L ) 

21 AAT2(L,K)=AAT2(K,L) 

13 CONTINUE 
15 II=JK*-1 
GO TO 12 
END 


00 

UJ 


OOOOOOOOOOOO 


SUBROUTINE DROP ( KO , NO V 1 , KE , NO V , NCS , FL/^G3 , FLAGG , P , XT , MTOC , KOFO , MUM ) 
COMMON X0(10) ,X1 (10) ,H(1 0, 10) , A(20, 10) ,GR0(10) ,BC(20) , 

IN 1 , NC , AAT2( 1 0 , 1 0 ) ,KOF , NC VO (20 ) , NLEC , NEC , MGSC , NGEX, KC 
DIMENSION NCV 1 ( 1 0 ) , NCV ( 1 0 ) , NCS ( 1 0 ) , P ( 1 0 ) , XT ( 1 0 ) , AMT (10), 

*UGR0(10) ,NTOC(10) 

REAL LAMDA 

INTEGER FLAGS, FLAGG, FLAG? 

FLAGGrI ALERTS CHK THAT A NSW POINT IS BEING RETURNED 
AND IT IS FEASIBLE SPACE. 

=2 ALERTS CHK THAT THE NEW POINT IS IN 
INFEASIBLE SPACE. 

=3 ALERTS CHK THAT THE NSW POINT IS IN INFEASIBLE SPACE 
AND ONLY THE EQUALITY CONSTRAINTS ARE KEPT 

FLAG7=1 WHEN THERE UNIT CONSTRAINTS PRESENT 

FLAG7=0 
KOFO=KF+KO 
B = 0 

DO 28 T =1 ,N1 
P(I )=GRO(I ) 

28 BrB+GROd )-^*2 
AGRO=SORT(E) 

DO ^3 T=1,N1 
43 UGRO(T )=GRO(I)/AGRO 
K = 0 

NT=NLEC+NEC+NGEC 
IF(KO.EQ.O)GO TO 3 
DO 2 T=1 , KO 

IF(HCV1 (I ).GT.NLEC. AND.NCV1 (I ) .LE. NLEC+NEC )FLAG7 = 1 
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NTOC(T ) = NCV1 (I ) 

2 CONTINUE 

3 CONTINUE 
IF(KE.EQ.O)GO TO 45 
DO 5 I=1,KF 

IF ( NC V (NCS ( I ) ) . GT . NLEC . AND . NC V (NCS ( I ) ) . LS . NLEC +N EC ) FLAG? = 1 
NTOC(I+KO)=NCV(NCS(I) ) 

5 CONTINUE 
45 CONTINUE 
C 

***** AMT IS USED TO SET THE SENSE OF THE NORMALS 

AMT>0 WHEN CONSTRAINT IS LE TYPE & GRAD POINTS TO 
FEASIBLE SPACE 

21 FLAG6=1 

HRITE(6, 330)MTOC(I ) ,1 = 1 ,N1 

330 FORMATCTHE CONSTRAINTS TO BE DROPPED ARE', 2013) 

AMVrl . 

IF(KO.EQ.O)GO TO 42 
DO 37 1=1, KO 

IF(NCV1(I).GT. NLEC. AND. NCV1(I).LE.(NLEC+NEC))GO TO 37 
B = 0. 

DO 36 J=1,N1 

36 B=B+A(NCV1(I ) , J ) »( X0( J )+UGR0( J ) * . 5E-5 ) 

C 314 FORMAT('I,NCV1(I),BC,B' ,2I3,2E13.6) 

TF(BC(NCV1(I))-B.LT.0.AND.NCV1(I).LE.NLEC)G0 TO 41 
if(bc(ncvi(i))~b.gt.o.and.ncvi(i).gt.nlec+nec)go to 41 

37 CONTINUE 
42 CONTINUE 

IF(KF.EQ.O)GO TO 4? 

DO 39 I = 1,KF 

IF ( NC V (NCS ( I ) ) . GT . NLEC . AND . NC V (NCS ( I ) ) . LE . NLEC +NEC )GO TO 39 
B = 0. 

DO 38 J=1 ,N1 

38 B=B+A(NCV(NCS(T )) , J)*(X0(J)+UGR0(J)».5E-5) 
IF(BC(NCV(NCS(I)))-B.LT.O.AND.NCV(NCS(T)).LE.NLEC)GO TO 41 
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00 

Os 


TF ( BC ( NC V ( NCS ( T ) ) ) ~B . GT . 0 . AND . NC V ( NCS ( T ) ) . GT . NLEC +NEC )G0 TO 1 

39 CONTINUE 
47 CONTINUE 

40 AMV=~1. 

FLAG6=2 

41 CONTINUE 

FROM 21 TO HERE IS USED TO SORT THE CONDITION 
WHERE THE GRAD POINTS INTO FEAS. SPACE FROM XO 
AMV=-1 WHEN GRAD POINTS INTO FEAS. SPACE. THIS IS 
KNOWN WHEN THE ONLY CONSTRAINT VIOLATED IS AN EQUALITY 

DO 26 1=1 ,N1 
P(I)=GRO(I) 

XT(I ) = XO(I ) 

DO 26 J=1,N1 
26 AAT2(I, J)=0.0 
K0F=0 

IF(FLAG7.NS. 1 )G0 TO 44 
WRITE (6, 318) 

318 FORMATCONLY THE EQUALITY CONSTRAINTS WILL BE KEPT') 

FLAG6=3 

RETURN 

THE REST OF THE SUBROUTINE CALC'S A VECTOR THAT IS THE 
AVERAGE OF ALL THE CONSTRAINT NORMALS AND IS GIVEN A 
DIRECTION OPPOSITE TO THE GRADIENT VECTOR. 

44 CONTINUE 
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SUBROUfWE PEASPT(PLA(x5) 

****-)k-**)^-llr*-)lr********)t**********i^***** ********** ********************* Mr*** 

Tills ROUTIifE BUILDS A SIMPLEX TABLEAU GALLED THE AUG MATRIX 
THE MATRIX IS CONSTRUCTED ROW BY ROW AND IS STRUCTURED AS FOLLOWS 
AUG = [ A : B ; C ] 

A - THIS IS TilE CONSTRAINT MATRIX WITHOUT X'S>0 
B - THIS IS A DIAG. MATRIX THAT CONSISTS Ox^’ THE POS. 

SLACK'S(NLEC) AND THE POS. ARTI?ICIALS(NEC ) . 

C - THIS IS A MATRIX OP ALh THE WEG. SLACKS ( .>fGEC ) . 

* iHfr *********** -iHHMt* ■)(• *»*•** ********************************* -X- **** * 

COMMON X0(10),X1 (10),H(10,10),A(20,10) ,GRO(1 0) ,BC(20) , 

1 N 1 , NC , AAT2 ( 1 0 , 1 0 ) , KOP , NC VO ( 20 ) , NLEC , NEC , N GEC , N GEX , KC 
DIMENSION C1 (30) ,C2(50) ,AUG(10,30) ,BCC(20) ,IPVTK(20),C:)(20) 

WRITE( 6 , 31 0 )NC , NGEX, NLEC , NEC , N 1 
310 PORMAT( 'NC, NGEX, nLEX, NEC, N1 ',513) 

INTEGER PLAG5,IPVTK 


***** xp THE NO. OP EQUALITY CONSTRAINTS 
DIMENSIONS. 


IS EQJAij TO The PROBLEM 


iP(NEC.i.T.N1 )G0 TO 2? 

22 WRITE(6,308) 

308 FORMAT (’TOO MANY EQUALITY CONSTRAINTS FOR 
PjjAG5 = 1 
RETURN 
27 CONTINUE 
NCI =NC-NGEX 
DO 4 1=1 ,NG1 
3CC(I)=3C(I) 

IPVTK(I)=0 


THE DEGREES 


OP FREEDOM ' ) 


LOAD A INTO AUG(A) ROW I 


DO 1 J = 1 ,N1 
1 AJG(I,J)=A(I,J) 


00 



oooooo ooo ooo 


00 

00 


LOAD 1 INTO AUG(B) IN THE PROPER DIAGONAL POSITION 

DO 2 J = 1 ,NC1 
AUG(I,N1+J)=0.0 

2 IFd.EO. J)AUG(I,N1+J) = 1 .0 
IF(MGEC.E0.0 )G0 TO 4 

LOAD -1 INTO AUG(C) FOR THE ARTIFICIALS OF GE CONSTRAINTS 

DO 3 J=1,N0EC 
AUG(I,NC1+N1+G)=0.0 

3 IFd.GT. (NLEC+NEC) .AND. (I-(NLSC+WEC)) .EQ. J)AUG(I,N1+NCUJ) = -1 .0 

4 CONTINUE 

DO 7 1=1 ,NC1+M1+NGEC 
Cl (I )=0.0 

THIS ALLOWS THE "CJ" ROW TO CONTAIN ONLY 0 EXCEPT 
WHERE THERE WAS AN = TYPE CONSTRAINT OR WHERE 
AN ARTIFICIAL HAS BEEN ADDED. THE "CJ" ROW IS 
TERMED Cl IN THE PROGRAM. C2 IS THE COST ROW. 

IF(I.LE.(N1+NLEC).0R.I.GT.(N1+NC1 ))G0 TO 7 
Cl (I ) = 1.0 
7 CONTINUE 
DO 30 1=1, NCI 
30 CB(I )=C1 (I+N1 ) 

17 CONTINUE 
WRITE(6, 301 ) 

301 FORMAT('AUG MATRIX') 

WRITE(6, 302)(N1 , ( AUGdI, JJ) , JJ = i ,N1 ) ,11 = 1 , NCI ) 

302 F0RMAT(NE10.3) 

WRITE(6, 302 )(HC1+NGEC, (AUGdI, JJ) ,JJ = 1+N1 ,NC1+MGEC+N1 ) ,11 = 1 , NCI ) 
C***** CALCULATE THE COST ROW 
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DO 6 T=1 ,NC1+N1+MGEC 
C2(I)r0.0 

J IS THE ELEMENT NO. IN COL. J 
DO 5 J=1,NC1 

5 C2(I )=C2(T )+CB(J)*AUG(J,T ) 

C2(I)=C1(T)-C2(I) 

8 CONTINUE 
WRITE (6, 330) 

330 FORMATCTHIS IS THE COST RO'W ) 

WRITE (6, 304 )C2 (LL ) , LL = 1 , NCI +N 1 +NGEC 
304 FORMATOOEIO. 3) 

***»* the COST ROW IS SCANNED TO DETERMINE TE IN FEASIBLE SPACE 

DO 9 1 = 1 ,NC1+NUNGEC 
IF(C2(I).GE.O)GO TO 9 

If A NEG COST IS FOUND THEN THE PT. NOT YET FEASIBLE 

GO TO 10 

9 CONTINUE 

«xxx« THIS IS A FEASIBLE PT . CANDIDATE 

GO TO 18 
10 CONTINUE 
1 = 1 

DO 11 J=1 , NC1+N1+NGEG 

THE COST ROW IS SCANNED FOR THE LOWEST VALUE 
THE T'TH COL. 13 THE PIVOT COL. 
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11 ■TF(C2(J ) .LT.C2(T ) )i:=J 

CALCULATE THE RATIO'S AND ETND THE PIVOT ELEMENT 
BY EJNDTEIG THE ROVJ WITH THE LOWEST RATIO 

BC12=1 .E19 
J = 1 

DO 12 L=1 , NCI 

TF( AUG(L,I ) .GT.O. )GO TO 2U 
BC1=1 .E20 
GO TO 23 

24 BC1=BCC(L)/AUG(L,T) 

2 3 CONTINUE 

IF(BC1.GE.BC12)GO TO 12 
J=L 

BC12=BC1 

12 CONTINUE 

IF(BC12.LE. 1 .E18)GO TO 28 
WRTTE( 6 , 315)1 

315 FORMATCCOLUMN' ,13, ' IS ALL NEGATIVE OR ZERO') 

FLAG5=1 
RETURN 
28 CONTINUE 
C 

IPVTK KEEPS TRACK OF THE ELEMENTS IN THE BASIS. IPVTK(J) = I 
C MEANS THE J'TH BASIS VARIABLE =THE I'TH VARIABLE 

C 

IPVTK(J ) = I 
C 

C»***Jt AUG(J,I) IS THE PIVOT ELEMENT AMD THE J'TH BASE VAR. 

C IS REPLACED BY THE I'TH VARIABLE 

Cfc***»f NOW NORMALIZE THE PIVOT ROW 
C 

BCC(J)=BCC(J)/AUG(J,I) 

AUGT=AUG(J,I) 
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DO 14 L=1 ,NC1+M1+NGEC 

14 AUG(J,L)=AUG(J,L)/AUGT 

the next tableau will now be DEETNED 

DO 16 L=1,NC1 

***** SKIP THE PIVOT ROW 

IF(L.EQ.J)30 TO 16 
AUGT=AUG(L,I ) 

BCC(L)=BCC(L)~BCC( J)*AUGT 
WRITE(6,314)L, J,BCC(L) ,BCC(J),AUGT 
DO 15 K=1 ,NC1+N1+NGEC 
AUG(L,K)=AUG(L,K)-AUG( J,K)*AUGT 
314 F0RMAT('BCC(L)=BCC(L)-BCC(J)*AUGT'/2I3, 3E13.6) 

15 CONTINUE 

16 CONTINUE 
WRITE(6,305)J,I 

305 FORMATC PIVOT ELEMENT IS ',21 3) 

AT THIS POINT THE BASIS VARIABLES MULTIPLIER COEFFICENTS 
ARE UPDATED. 

C3(J)=C1(I) 

GO TO 17 
18 CONTINUE 


***«* THIS IS A FEASIBLE PT. CANDIDATE 

***** ALL BCC'S WILL BE SCANNED TO BE SURE THEY ARE POSITIVE 

DO 19 1 = 1 , NCI 
IF(BCC(I).GE.O.O)GO TO 19 
GO TO 21 
19 CONTINUE 


'O 



VO 

K> 


c 

SET THE INITIAL FEASIBLE POINT 
C 

DO 20 I=1,NC1 

20 X0(I)=0.0 

DO 29 1=1 ,NC1 

IF(IPVTK(I).GT.0.AND.IPVTK(I).LE.N1 )X0(IPVTK(T))=BCC(I) 
29 CONTINUE 

WRITE (6, 320) 

320 FORMATCTHIS THE INITIAL POINT') 

WRITE (6, 321 )N1 ,(X0(I),T = 1,N1 ) 

321 FORMAT(NE13. 6 ) 

RETURN 

21 WRITE (6, 300) 

300 FORMATCHO INITIAL FEASIBLE POINT FOUND') 

FLAG5=1 

RETURN 

END 
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